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A digital computer program, originally developed by George E. Fink, 
capable of determining the in-and/or out-of-plane vibration frequencies 
of a single plane piping system, using the basic transfer method, is 
modified and augmented by the writer. The program is modified to use the 
writer's method to reduce the amount of calculation and is augmented to 
use Marguerre and Uhrig's modifying frequency determinant method and 
Pestel and Mahrenholtz ' s remainder method for higher natural frequencies. 

An analysis is made of difficulties encountered with the original 
transfer method and utilizing the modified methods. 

The program accepts branched systems and non- stiff intermediate 
supports. A discussion and an explanation of how to use the program are 



included. 
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CHAPTER I 



INTRODUCTION 



1-1 General remarks. 

Today, in increasing measure, natural frequencies of a system capa- 
ble of linear vibration are being determined with the help of the method 
of transfer matrices. This method is well explained in the recent book 
by E. Pestel and Leckie [3].* 

This method which has been known for more than ten years, proceeds 
with particular simplicity especially with the use of digital computer. 
However, it meets with important numerical difficulties. 

Numerical difficulties were observed by George E. Fink in his U. S. 
Naval Postgraduate School Master’s thesis ’Vibration analysis of piping 
systems” [4] . He used the straight-forward method of transfer matrices 
for calculation of the natural frequencies of vibrating piping systems. 

The purpose of this paper is to deal with numerical difficulties 
and with the remedies which have been developed and also to investigate 
reduction in calculation. 

1-2 Scope of work. 

Numerical difficulties are encountered at higher natural frequencies 
and for systems having very stiff exterior springs. In examining these 
difficulties the paper by Marguerre and Uhrig was encountered [2]. This 
tends to explain difficulties both with the original transfer method and 
with the Pestel and Mahrenholtz's modified methods [1]. Although, 
Marguerre and Uhrig suggest several different methods for dealing with 

★Numbers in brackets refer to the bibliography, page 37. 
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the difficulties, those which appear to be certain to eliminate the 
difficulties are so far outside the spirit of the transfer method that 
the advantage of this method would be lost. Therefore, it was decided 
to proceed with the transfer method to get as much good from it as one 
could. 

The writer set out to investigate the remedies for failure at high 
frequency. The digital computer program "VIPIPE" presented in the Ap- 
pendices was used to determine the natural frequencies of in-and out-of- 
plane vibration of a planar piping system based on the transfer method. 

This program was originally developed by Fink [4], but has been modified 
and augmented by the writer. 

1-3 Assumptions and limitations. 

1. System is conservative. 

2. Distributed mass or/and lumped mass model is used. For the lump- * 
ed mass model used in curved section (or in straight section) sufficient 
subsectioning is carried out to give acceptable answer. 

3. The transfer method surely fails if the system has rigid exterior 
elements. The hangers are introduced as a support by a point matrix with 
appreciable linear and/or torsional compliance. 

4. Program VIPIPE is limited to a system of a maximum of fifty com- 
ponents, twenty branches, and twenty hangers. It can only determine in- 
and/or out-of-plane natural frequencies of planar system. 

5. Boundary conditions should be such that half of the state quan- 
tities are vanishing.* 

*Refer to Appendix B. 
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1-4 Notation 



[ ] 
{ } 


matrix. 

column matrix. 


X,Y,Z 


cartesian right handed coordinate system. 


u,v,w 


displacement in the XYZ direction respectively. 




rotation about XYZ axis respectively. 


M 


bending moment. 


N 


normal force. 


T 


torque. 


V 


shear force. 


Z. 

1 


state vector at location i. 


'V 


state matrix at location i. 


l u i«! 


transfer matrix from i to i + 1. 


A(U» 


frequency determinant. 




modified determinant.* 


[N] 


non-vanishing state matrix of Z . 

o 




element of[U]. 


[S] 


spring matrix. 


[S] 


purged frequency determinant in the P-M method.* 


[I] 


unit matrix. 


ro] 


null matrix. 


d 


displacement vector. 


[d] 


displacement matrix. 


P 


internal force vector. 


[p] 


internal force matrix. 



*In the text, the following abbreviations are used. M-D method stands 
for Modified Delta method and P-M method for Pestel and Mahrenholtz' s 
remainder method. These will be explained in what follows. 
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(GJ 


coordinate transfer matrix. 


M, M, (R3j 


square sub-matrices . 


v 0 


column vector of every elements value of 1. (tf* p* 


R(oa) 


remainder. 


X 


correction factor vector. 


Po > P 1 > Pi 


assumed state quantities at the starting boundary. 


Xq , X, , x a 
Y 0 , Y, 


correction factors in the P-M method. 


D| , Dj 


purging factors in the M-D method. 


A. 


eigenvalue. 


X 


eigenvector. 
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CHAPTER II 



TRANSFER METHODS AND NUMERICAL DIFFICULTIES 



2-1 The state vector and transfer matrix. 

A. State vector. 

The state of vibration at location i, specified by state quan- 
tities, is described by the state vector Z^. This vector is a column 
matrix whose elements completely describe the instaneous displacements, 
both rectilinear and angular, of that point from its quiescent position, 
as well as the corresponding rectilinear and angular forces in the member 
at the same point and at the same instant, i.e., those forces, if a cut- 
ting plane were passed through the point at the instant of interest, would 
have to be applied to each cut face to prevent relative motion between 
them. 



In the case of planar piping system, the state vector has six ele- 
ments, three pertaining to displacement, and three pertaining to force. 

The state vectors, for the in- and out-of-plane cases respectively, 
as used in VIPIPE, are: 

u displacement in X direction, 

v displacement in Y direction. 

U slope in X-Y plane with respect to X asix. 

Mj moment about Z axis. (2-1-1) 

V|. shear in Y direction. 

N normal force in X direction. 



z 4> " 
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twist about X axis. 



Z of 



( 2 - 1 - 2 ) 



s 

T torque about X axis. 

-w (negative of) deflection in Z direction. 

•jJ slope in X-Y plane about X axis. 

My. moment about Y axis. 

* shear in Z direction. 

' o 

Thus, the state vector in this application is of order 6x1, while the 
system state matrix* is 6x3 and the transfer matrix of a section is 6x6; 
however, in developing theory we will use the more general case 2rxl, 2rxr, 
and 2rx2r respectively. A local XYZ right handed cartesian coordinate 
system and the elements of the state vector are arranged as in Fig. 2-1-1. 




Fig. 2-1-1 

B. Transfer matrix and transfer method. 

From the state vector at location i (node i) , the state 



*The concept of "state matrix" is introduced in Section 2-3. 
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vector at location i + 1 is given by 



z i + l " K+J * Zi. (i “ 0,1,2,..., n) (2-1-4) 

Here, matrix relating the state vectors of node i and node i + 1 

is termed a transfer matrix. 

We number the node i + 1 immediately following the node i from the 
left end, the beginning of the main system. From the recurrence relation 
2-1-4, one easily obtains 

z „ ■ ( u nJ W- • • N • (”2) • M • Z o 

- (u) . Z o (2-1-5) 

which is a linear relationship between the state quantities at location 
0 and n of the system (boundaries of the main system). We call this 
method the transfer method. 

Since of the 2r state quantities at each boundary r quantities are 
vanishing, only r are established by the boundary condition at location 
0, and r homogeneous equations are established by the boundary condition 
at location n. The natural frequencies of the system are given by the 
zeros of the determinant of the corresponding r homogeneous equations. 
Since the application of interest is to the determination of the natural 
frequencies of piping systems, it would be instructive to illustrate the 
origin and nature of frequency dependence in transfer matrices. For thjs 
end, let us construct a point matrix for the in-plane case connecting Z i 
with (z£ and Z^ are state vectors of right side and left side, 
respectively, of a point mass m^ at node i) . 

To begin we assume that all of the elements of the state vectors 
vary with time, each having the factor cos . 



11 



Noting that the deflection, slope, and moment* are continuous across 



the concentrated mass m^, so that 




From the free body diagram of Fig. 2-1-2A, B, the following relations 
are established. 



c 

rr 

N-/ 

II 

C 

• 

o 

o 

CO 

ft 


— * -U u) coslO t 


(2-1-6) 


< 

/-N 

rr 

V-/ 

II 

< 

O 

O 

CO 

& 

ft 


= -vu) cosOt 


(2-1-7) 


R \ u 

V i • cos W t-V^_ cos 




(2-1-8) 


N ^ • cos U) t-N^ • cos 


U) t - - 7T1- 


(2-1-9) 


and substituting Eqn. 2-1-6 into Eqn. 2-1-8 and Eqn. 2-1-7 


into Eqn. 


2-1-9 we have 






> 

H 

> 


+ mv uf 


(2-1-10) 


H? -HV 


+ mu i^ 2 " 


(2-1-11) 



Now, the relation of the state vectors of the right and left side of 
the point mass m in matrix form is: 

2 i = [ u p] ‘ z i 

i. • 6 • ^ 



/ > 
u 


= 


1 


0 


0 


0 


0 


0 


• 


u 


V 




0 


1 


0 


0 


0 


0 




V 






0 


0 


1 


0 


0 


0 




v* 






0 


0 


0 


1 


0 


0 










0 


m ^ 


0 


0 


1 


0 






N 


LR 


mC) i 


0 


0 


0 


0 


1 




N 


S J 


<* 














/ 



[UpJ is a point matrix for in-plane vibration. It neglects any rotational 
inertia that the point mass may have; however, rotational inertia could 

*The moment is continuous only if ehere are no point elements having finite 
rotational inertia. 
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be accounted for by including anothoranonfceto* term in the matrix. 




> f 







inuf u 



* 






-rnvOV 

Fig. 2-1-2A Fig. 2-1-2B 

The derivation of transfer matrices is treated in great detail in refer- 
ence [3] and transfer matrices for various components of planar piping 
system are incorporated in VIPIPE. 



2-2 Failure of transfer method. 

For higher natural frequencies and for systems having very stiff 
exterior springs numerical difficulties are encountered. In these cases, 
the numerical value of the frequency determinant is given as the differ- 
ence between large numbers. So practically the natural frequencies can 
no longer be determined. 

From a mathematical point of view, in case the transfer matrix [U) 

A? 

has a single dominant real eigenvalue, the columns of [U] become more 
and more parallel as k increases, and approach the first eigenvector x^ 
[6] [ 2 ] . In this case the transfer matrix degenerates. 

[Uj\, - K> x, < 2 - 2 ' V > 

This could be the case for the lumped mass system as the number of lumped 
masses increases. This was observed in the straight section lumped mass 
system, but for the curved lumped mass system this was not observed. 
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As the frequency increases (elements of some transfer matrices are 
functions of frequency), the same phenomenon was observed in the distri- 
buted mass straight section transfer matrix. For example, for the sample 
problem of Appendix E-2-1, 

At frequency of 20 rad/sec, 

2 

the largest eigenvalue: = .11484022x10 

second largest eigenvalue: \ = .99980414 + jO. 01979160 

transfer matrix (U) is: 



.99980413 


0 


0 


0 


0 


.90535567V" 5 * 


0 


.25105587V 1 


. 25986435^ 


.22549093V" 2 


. 14255229V 


0 


0 


. 3084007 3\ _1 


. 25105587V 1 


. 26651055 V4 


.22549093V" 2 


0 


0 


.47565399\ 4 


.30070198 x6 


•25105587V 1 


. 25986435V 3 


0 


0 


.56219618V 


.47565399V 4 


.30840073 V" 1 


. 25105589V 1 


0 


-.43265609\ 2 


0 


0 


0 




.99980413 



*V 5 means 10 5 , etc. 

In this case, that is, at this frequency, we can see that there is no 



phenomenon of the columns becoming parallel. As the frequency Increases 
this single dominant real eigenvalue keeps increasing, and at frequency 
of 1700 rad/sec. 

the largest eigenvalue: \, = .5937220 x lO 1 ^ 

second largest eigenvalue: \ 2 , = -.11136847 + j. 99377918 

transfer matrix [u] is: 
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-.11136847 


0 

. 14843074\ 10 


0 


0 


0 


.53482189 v5 


0 


.13191200 \ n 


.12023306\ 5 


.10685242 x6 


0 


0 


.16701804\ 9 


. 14843074 \ 10 


.13528929\ 4 


.12023306'''** 


0 


0 


.18324149\ 15 


.16284870\ 16 


. 14843074\ 10 


.13191200 N ' 11 


0 


0 


. 20618799\ 14 


.18324 14 9\ 15 


. 16701804^ 


. 14843074\ 10 


0 


-.18465906\ 6 


0 


0 


0 


0 


.11136847 



Normalized values of each column of [u) are: 



.603103^" 7 


0 


0 


0 


0 


-. 4802273^ ” 5 


0 


.8100280\‘ 6 


.8100279 V 6 


.8100280\' 6 


.8100280 V 6 


0 


0 


. 911464 l \' 7 


.91146407 V 7 


. 911464 lV 7 


.91 14641 V 7 


0 


0 


1 


1 


1 


1 


0 


0 


.1125225 V 1 


.1125225 V 1 


.1125225\ _1 


.1125225 V 1 


0 


1 


0 


0 


0 


0 


1 



We can see several of the columns are almost parallel to each other. From 
this, it is easy to see the numerical difficulties we should encounter. 

This depends on the nature of the transfer matrix. For systems made of 
various sections, multiplications of different transfer matrices may result 
in the difficulties taking care of themselves. 

2-3 Reduction in calculation in matrix multiplication. 

In equation 2-1-5, [u] = [Uj . . . . [Uj . [U,] > so we 

should perform matrix multiplication of 2rx2r matrix times 2rx2r matrix 
successively. 

Now using the fact of the boundary condition, introduce Z Q as 
2rxr matrix called the "state matrix." This is formed from the state 
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vector as follows, all elements of the state matrix are zero except that 
in its ith column, unity appears in the same position as the ith general 
ly nonzero term reading down the state vector. Then we can multiply 
2rx2r matrices by 2rxr matrices and thus can reduce the calculation 
considerably. 

For example, in the in-plane case, fora built-in left end 



Z 



o 



' 0 0 O' 

0 0 0 

0 0 0 

1 0 0 

0 10 
0 0 1 , 



The frequency determinant is obtained in the usual way from the boundary 
conditions at the right end. 

We use the term "Delta method 11 to refer to this method which was 
discovered by the writer during the course of this study. This pro- 
cedure is used in what follows instead of the usual procedure of the 
transfer method. 



2-4 Pastel and Mahrenholtz 1 s remainder method. 

Instead of the frequency determinant, Pestel and Mahrenholtz intro- 
duced a remainder [l] # [3] . This method, presenting the system state 
matrix as 2rxr matrix, reduces the amount of calculation considerably 
compared to the original transfer method, and claims greater accuracy at 
higher frequencies because the remainder, which is. the difference between 
small numbers, is itself small. This will be explained later. 

A. Remainder. 

For the non-vanishing state quantities of the state vector Z q 
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which depends on frequency \0 we use estimated values which approximate 



the true values, following which the (i = 1, 2, 3, ...,n) are deter- 
mined. 

In the state vector Z 0 , since for eigenvibration only the ratios of 
state quantities are of interest, one of the non-vanishing state quanti- 
ties is assigned the value of 1. 

For example, for the in-plane case, at certain frequency, a canti- 
levered pipe which is free at the left end and fixed at the right end has: 

.o 



7 

**0 



O 0 

p. + x; 



0 O 

p* + x * 



(2-4-1*) 



0 
0 

^ 0 

(in what follows, superscripts 0 and 1 denote the round of calculation; 
i.e., super-script 0 means first round of calculation with approximated 
state quantities in Z ft , and superscript 1 means second round of calcula- 
tion with improved state quantities in Z fl ). In equation 2-4-1', we 

assigned the value of 1 to the non-vanishing state quantity in the first 

o o 

row of Z 0 , and the quantities of P, , P z are the previously chosen 
approximations and X* , X° are correction factors which are unknowns. 

The state vector Z° of (2-4-1') in other matrix form is: 






0 

0 

o 
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0 

1 

0 

0 

0 

0 



0 
0 
1 
0 
0 
0 J 



1 

*: 

x! 



(2-4-1") 



For the actual matrix multiplication Z 0 of (2-4-1' 1 ) is used. After a 
round of calculation, the approximated state quantities P° , P° will 
be improved in the form 



p; = p? + 

P 2 = K + K (2-4-2) 

where X° and "X ° are corrections obtained by this process; the 
difference between and X^ , and will be shown later. 

Then, the improved state vector z' 0 will be: 




0 

0 

. 0 - 



When we start first, usually P° , P° 
absence of better information. From the 
Zq of (2-4-1"). 



(2-4-3) 



are set equal to 1, in the 
same method as (2-1-5) using 



Z?I = [Ut.] • (U-n-,] • (u,j • z 0 ° 

From the boundary condition at location n_, i,#. , from the three vanishing 
elements of , we get 

0 * A(U>)’ (N°) • x°= (S]‘ X° (2-4-4) 

[N] : nonvanishing sub-matrix of [Z 0 ]. 

X : correction factor vector. 

A(<a)) : frequency determinant. 

[S] : A(u>)-(N) 



For the cantilevered pipe of the example in the in-plane case, the 
frequency determinant A(u)) , see Ref. [3] , consists of the top left 
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3x3 submatrix of the transfer matrix [U] . 
Equation 2-4-4 may be written: 



'O' 


at 


'u„ 


U,1 


u.13' 


• 


1 


0 


N 

0 


• 


1 ■ 


0 




U2I 


U 2 2 


Uj 3 




Tf 


1 


0 




x; 


0 

k / 




U 3 I 


U32 


%3 , 




°Q? 

/ 


0 


1 ( 




V. 

v ■ — 



SB 


’U,| 


+ ua p; + Up Pi 


U|2 


(Jl|3 




Urt 


+ + UnPl 


lAii. 






. u 3t 


+ u j 2 p; t U 33 p; 


U 32 


U 3 3 < 



• 


t \ 

1 


5 


X 

— 0 


/ 


P;J 



as 


' 5 11 


5|Z 


' 


• 


f * * 

1 












x; 




$31 

k 


632 


5 S3 

/ 




[xlj 


4-5, the 


first 


column 


of [j 



(2-4-5) 



numbers, but the second and third columns remain unchanged. We call 
this "purging" the first column. From equation 2-4-5 we may write, 



V 


— 






% 








5^2 


^3 




x; 

s / 


.^31. 




v %2 


$33 

y 







(2-4-6) 



Equation 2-4-6 has first, three equations with two unknowns (two cor- 
rection factors) , and second it has the purged column to the left and 
unpurged columns as coefficient matrix. 

Equation 2-4-6 is the same as; 



U.|| 


= 


Ul2 


U| 3 


P21 






U 2 3 


/—r— 

£- 

<Al 

\ 




u 32 


£■ 

Uo 

N 



p:+ x; 

p;+ x; 



(2-4-7) 



From 2-4-7 we choose any two equations. If we choose the first two, 
then by Cramer's rule; 
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— 


U11 


U13 




— 


Uu 


U.11 






Uai 




p° _ A* _p° 


x;-- 




Ua| 


R. - A) -^<2-4-8) 


X,- 


U |2 


Vi 13 


A, *■ 


U| 2 . 


U |3 






**.23 








U23 





The A 's are defined by the preceeding expressions. If U)=U)j^, 
where is a natural frequency of the system, then the third equa- 

O 0 

tion will be satisfied if we substitute these values for X, and . 
Otherwise, substituting the above value for X° , we get, from the 
third equation, a somewhat different value for X° ; we denote this by 
the symbol Y° 



Y 



o 

X 



1 * 3.2 Ax 1 * 3 1 

U 33 At U33 




(2-4-9) 



We define the remainder R°(lO) to be the difference between these 
values. 



R # (*)= x;- Y 



o 

2 



A3 U33 U31 

A I U33 A| U35 



(2-4-10) 



As we see from equation 2-4-10, for fixed 0 , the remainder is indepen- 

0 O 

dent of the assumed state quantities and P a . If the remainder is 
zero, ^ = u)j^ since all the prescribed boundary conditions will be 
satisfied simultaneously. 

For the next round of calculation for the same frequency, form the 
arithmatic average of the two values X and Y for the correction factor 
% , i.e. , 

= X°, 

Then, new starting values are: 
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A2 

A\ 




1 ( a 3 


I U 32 Aa 




2 l A, 


VU 33 A| 


U3 iJj 



(2-4-11) 



O O 

which are also independent of P, and P 2 . A repeated round of cal- 
culation gives 

x;= x! = -p;- ~ = o 

-r « - Yi) = 0 

R'(0)= xj - Y* = RVx> 

therefore, 

x; = -Yi - 4- Rio) 

Now, we see that the iteration is completed after the first round of 
calculation and the remainder is no longer the difference between very 
large numbers, because the remainder is of the same order of magnitude 

i i 

as the corrections X 2 and Y 2 . During the iteration seeking the natural 
frequencies, we calculate R°(^i) and pj_ ( tOj ) 's (i=l,2), and increase 
frequency by an amount AO to get lOa-0|+AlO. Then get R°(lOj) using 
pj( v) 2 ) = p] (i*)|) (i = 1,2). The above example is for the case r = 3. In 
other cases for r ^>1, the spirit is the same (in case r « 1, we cannot 
use the P-M method) . 

B. Various ways to get remainder. 

There are many ways to calculate a remainder. Each one has a 
different value at the same frequency (ti)^u)jg) , but passes through zero 
at O — . 

Some fail and some give good results (see section 2-4-D) . 

The possible ways arise from the following choices; 
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1. Which one of the nonvanishing state quantities will be 
given the value of 1 in the state vector Z^. 

2. Among the r homogeneous equations, as in equation 2-4-5, 

o o O 

which r -1 equations will be chosen to get X, * , ...X r _ j . 



i 

3. Which of the X° 's (i 1,2,..., r-1) in the rth equation is 
to be regarded as Y° and solved for in terms of the other X° 's‘. 

We will call these different ways by the name "sub-procedures." 



C. The relation between a remainder and the frequency determinant. 
From the equation 2-4-10, the remainder is 

r( ^ ■ ^ 33 ^ 34 - U32 42 lb|A| 

U33A1 



From equation 2-4-8, we see that A's are the minors of an element of 
the frequency determinant 



A|=M(W3l) , = M( “32) , A3=M(U 33 ) 

where M(u 3 j ) means the minor of 

The remainder is 

A(l0) 



R< 0 ) = 



U33 A 1 



where A( 0 ) is seen to be the frequency determinant. 

Expressing a remainder in general in terms of the frequency deter- 
minant; 



R( l*) ) 



A(0) 

Utf Ars 



(2-4-12) 



Uj,^. : coefficient element corresponding to Y ° of rth row. 

A (-5 : minor of u hs in A ( 0 ). 

Up. g : element in rth row and the column corresponding to state 
quantity of 1 in Z Q . 
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D. Failure of the P-M method. 

In case the denominator of the remainder is zero, the P-M method 
certainly fails.* (Refer to equation 2-4-10 and 2-4-12). When the deter- 
minant in the denominator becomes nearly equal to zero or changes sign 
passing through zero, the remainder is difficult to evaluate: at certain 

frequencies the remainder may show an infinite or jumping phenomenon. 

However, this does not necessarily show up at the same frequency 
with different sub -procedures, because this phenomenon depends solely 
on the nature of the transfer matrices. 

Figure 2-4-1 and Figure 2-4-2 show that while one sub-procedure 
shows an inifinite phenomenon at certain frequency, another sub -pro- 
cedure does not show the same phenomenon at this frequency. 

Pestel and Mahrenholtz illustrated their method with the case of a 
homogeneous beam (2r =4). It leads to very good results for the bend- 
ing frequencies. The reason is, because r = 2, there is no phenomenon 
of indefiniteness causing failure mentioned above and due to purging 
effect. The P-M method does show improvement also for the case of r great- 
er than two. 

For sufficiently high frequencies, the P-M method fails. The reason 
is as follows, if the columns have already become parallel, purging does 
not do any good, and if round-off error reaches critical point, where 
remainder is no longer reliable, this method fails. 



*When this is the case, program VIPIPE skips to the next sub-procedure. 
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Figure 2-4-1 shows fourth and fifth mode frequencies of System 5* at 
frequency 826.88281950 rad/sec, and 1126.83799669 tad/sec (P-M method- 
subprocedure A with double precision arithmetic). Figure 2-4-2 shows 
infinite phenomenon at frequency 939.00298809 rad/sec while Figure 2-4-1 
does not show this phenomenon at this Frequency. And Figure 2-4-2 has 
print out-put with mode frequency 826 88281915, also with 939.00298809 
rad/sec (P-M method- subprocedure B with double precision arithmetic). 

From Figure 2-4-2, we see that the latter is a false mode frequency and 
the former does not show upon the graph due to the infinite phenomenon 
of the latter, but we know that there is a mode frequency at 826.88281915 
rad/sec (this also can be checked with the frequency and remainder value 
print output) . 

2-5 Modified Delta method. 

Instead of purging just one column vector among r columns (as in the 
P-M method) of frequency determinant, Marguerre and Uhrig showed in their 
paper [2] a possible way of extending the purging concept for problems 
greater than 4th order. 

Through exactly the same method as Delta, we get a frequency deter- 
minant. At high frequency, because the columns approach parallelism, 
the value of the frequency determinant tends toward zero. If the columns 
are exactly parallel to each other, then of course all methods fail. If 
they become nearly parallel, then with a finite capacity of handling 
numbers, the results become meaningless. 

Using Maguerre and Uhrig' s purging concept, with the help of one 
column of the frequency determinant, obtain r-1 purging factors to purge 

*See Appendix E-3. 
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another r-1 columns.* For the case of sixth order (r*3) as used in 
program VIPIPE, the frequency determinant at a certain frequency is: 

A(ti>) =• Un U.|2 Ub 

^21 U22 ^23 

U 3 | Lhg U33 

Purging factors got with the help of the third column of A(lO) are 



D| : 



J_ ( _LkiL + 

_3 ^ Uo U33 ) 

I ( ^13- , Mjq U32 \ 

a V u ,3 T U *3 U33 ' 



(2-5-1) 



Modifying the non-vanishing sub-matrix of the starting state matrix [Z 0 ] , 
we have 



[N] = 



1 

0 



0 

1 



; D i -°2 



0 

0 

1 



(2-5-2) 



— 


U|| 


D,U |3 


Ua 


Dj U13 


U |3 




Uit 


D|U« 


U22 




u 23 




u* 


DiUi3 


U32 


Pj ^33 


U33 



Then, for the same frequncy, repeat a round of calculation to get the 
modified frequency determinant. 



(2-5-3) 



From the properties of determinants, 4 (^ ) is the same as (£>). 

But the first and second columns of 4^ (0 ) are not the same as those 
of A(u)) any more. If the columns of the original frequency deter- 
minant were exactly parallel, the modified columns (first and second 
columns) of the modified frequency determinant will turn out to be all 



*We use the term "Modified Delta method" to refer to this method. 
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zeros. In case the columns were nearly parallel, the first and second 
columns will turn out to be very small numbers, i.e., we will have purged 
two columns. So we can expect a more accurate value of the frequency 
determinant and can go further to higher natural frequencies. 

If = D 2 = 0, it degenerates to original frequency determinant, 
so it does not do any good. If = E , it likely means that first and 
second columns are the same. If this is the case, this method just fails. 

In some cases, especially systems of many sections, the frequency 
determinant does not show parallelism. Experience shows that this method 
gives good results in most problems. However, this method also fails 
at sufficiently high frequencies. 

2-6 Intermediate conditions - Branched system. 

The effect of a joint point can be introduced into the calculation 
by means of a single point matrix. 

This section is mainly for the derivation of the branch point matrix 
of the Modified Delta method and the Pestel-Mahrenholtz method. 

The derivation of branch point matrix for the straight forward ^trans- 
fer method (and thus also for the Delta method) is fully explained in 
Ref. [3] and [4] , but a brief discussion will be given here for the 
purpose of refreshing the reader's knowledge and moreover explaining the 
procedure for the M-D method and the P-M method. 

A. Delta method. 

Discontinuities are introduced in the forces whereas the main 
member deflection on each side of joint are the same, and identical to 
those of the branch at the same joint. 
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By a suitable coordinate transformation, it is possible to express the 
force system of the branch in terms of the main system coordinates. 

A free body diagram for the forces at the joint point is shown 
below.* 



1 JL 








Fig. 2-2-2B 



The forces at the joint in matrix form are. 



M' 
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’ M ' 
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'M 


V 




V 
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N / 




M , 


JL 


. M . 



TIB 



M 



TR 






■7R 



( 2 - 6 - 1 ) 



The displacement is continuous at the joint, so Z can be obtained by 
simply adding the force components of branch state vector to Z ^ . 

The column matrix of nonzero component at the starting boundary will be 
symbolized by V Q , Here V 0 = ^1, 1,1 J 



*Symbol /\ indicates that the vector elements are expressed in terms of 
the branch coordinate system. 
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(2-6-2) 
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From equation 2-6-2 

A "I - 

P-nB = U^. 2 )* t^i] * (2-6-3) 

The matrix product j ^ is called a spring matrix. To transform 

the coordinate system of the branch into main member coordinates at the 
joint, we use 









(frlj * d T)B 






P H6 ~ 


• P-rc 


here; 


G l = 1 
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(2-6-4) 
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where <p is a branch joining angle; refer to Appendix B-5. 
then; 



P,„- <U 

— ( 5 ) • d.n B 



(2-6-5) 



((») = (■%)• ( R 2)-( R l)'-[°l!! 



(sj is also called a spring 
matrix. ) 



Recalling the relation at the joint 

d TR ’ d JL “ d TIB 

p _ ■ p +P_o = P + r S 1 . d R 

‘tl r 71 B *JL 1 J 

we get in matrix form 



29 



■J*. 



' 3T_ 



i O'! • f d 

,s ij [p 

or = tU p ] • 

Here [Up] denotes the branch transfer point matrix. 



( 2 - 6 - 6 ) 



B. Modified Delta method. 

(In this section, the use of an asterisk * will denote a 
purged matrix or vector. Also note that the quantities d and p, with or 
without * and , represent 3x3 matrices whereas in the preceding sub- 
section they represented 3x1 vectors.) 

After we get the purging factors from^(u>), the new starting matrix 
at location 0 is 

[zh “ 1 0 ] , l H"J - 



N 



1 

0 

-D. 



0 

1 

-D, 



At the joint 



(O 
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d* ' 


f Z TU*] » 
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.P*. 




P . 



0 
0 
1 

(N°J 



TL 



For a branch, let the non-vanishing state matrix of the new starting 

state matrix [Z 0B ] be 

[N e ] = [100 

0 10 



-D, b -d 
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Here D, e , D^g are branch purging factors, the calculation of which 
will be explained later. At the joint 
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[R^J , [I^Jf [R*^] & fR* 2 ^ are defined by this expression. 

Just as in equation 2-6-5, 

[«*) - (&»)•( r/HrT]'' (&,) 

= [Gil - (R j)' (Mb) ' (Ns) ’ ( Ril'tfriJ 

* [&.]• (Rj- [ Ri] [&,) 

• (Si 

The spring matrix turns out to be the same as in the Delta method. 



At the joint 

p? e - (pj-(M-)-v.+ (p;j-v. 



- ( p.T ) - v 0 + (sHcUJ-v. 

[Sj is independent of(N g ] , so p* B is independent of [N 0 ] 



For Z^* we have 



Z TR* 



I 

s 



0 

0 



Z l = f U pl* Z ?u 



We see that branch transfer matrix [u p J is the same as in equation 2-6-6. 
Purging factors can be obtained from [R^] or The phenomenon of 

columns approaching parallel could occur in [Rj and [R^] ; VIPIPE com- 
putes purging factors from (R^j , because the writer was inclined to 
think that getting purging factors from [R^J would do more good in in- 
version of [Rj] . However, [R^J could have been used. 

Purging factors D (g and D 3 g are obtained in the same way as in 
equation 2-5-1. 

C. Pestel-Mahrenholtz method. 

At the joint the equilibrium condition must be satisfied, 
and the deflections should be continuous. 
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Here, [N* B ] is the nonvanishing state matrix of [ L^ e>j and 
branch correction vector, 
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Then, the state vector at locatior OB for the built-in, in-plane case 
is 



At the joint 



Z° = 

OB 



0 

0 

0 

p oB x; 

D° Y° 
i IB AqB 



Pi X 



0 

oe> 



mb 






z JL = 



A *■ 
^•TlB ~ 



• [NJ°J‘X C 



JU 



•(Ny*xi = 






x c 



(2-6-7) 
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where [Rj*J and^R^*] are 

(R?)= [R.H< e ) [Rll- [RJ-fN^J 

For[N°) and(X°J , refer to section 2-4. Also [R^J , [R^j are expressed in 
the same way as the preceding sub-section B. 

(S*J - (a-i)- [ rD- [R n^Cfr.l 

■ IGjMr* ) ‘(m*]- (nJJ-(r J • (<fi) 

- [6iHM-(R. f'(&i) 

- is) 



The spring matrix also turns out to be the same as in the Delta method. 



From the requirement for continuous deflections at the joint we get 
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(2-6-8) 



and 
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For Z JR we have 



Z 



7R 



t Sj * d - (S) • [d T J-(N°J • X° 

(p£)-x° + (s H d £j x ° 




£°p) Z JL 



We see that the branch transfer matrix [U p ] is also the same as in 
equation 2-6-6. From equation 2-6-8, the branch correction vector is 

X ob = ( R ?H G lH d £j*X° (2-6-9) 

After we get the main member correction vector (X ) , we can get branch 
state correction factors in terms of X° values. As explained in section 
2-4, iteration is performed in the same way. 
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2-7 Computer. 

In FORTRAN 63, floating point quantities have an exponent and frac- 
tional part. In single precision arithmetic, the fractional part has 
36 bits and the precision is approximately 10 decimal digits. In double 
precision arithmetic, the fractional part has 85 bits and the precision 

is approximately 25 decimal digits. In both cases the range of numbers 

, , ,.-308 ..308 

is from 10 to 10 

In the FORTRAN 63 version of VIPIPE (as described later), either 
single precision (for conservation of computer time) or double precision 
(for greater computational accuracy) can be selected. In the FORTRAN 60 
version, only single precision is available. 

The natural frequencies of a planar piping system are found using 
program VIPIPE. Program VIPIPE has been tested with FORTRAN 60 and 63 
using digital computer CDC 1604. The number of natural frequencies that 
may be found using VIPIPE and the accuracy, especially at high frequency, 
depends upon the floating point significant figure capacity. And the 
significant figure capacity required for any system is a function of 
the characteristics of the system itself. 
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CHAPTER III 



DISCUSSION 

3-1 Accuracy of methods and assurance of the solution. 

Confidence in accuracy of methods was established by comparing 
frequencies of several systems for which these frequencies could be 
obtained by other means, i.e., by means of either a closed form analyti- 
cal solution or by recourse to the Ref. [5] . 

The accuracy is very good for the lower natural frequencies in the 
values got using the Delta, method, the Modified Delta method, and the 
Pestel-Mahrenholtz' s method. 

In most cases, up to fourth mode frequencies, the difference be- 
tween the system natural frequencies obtained by analytical means and 
those obtained using VIPIPE are less than 0.004%. As frequency in- 
creases we got more accurate natural frequency values from the M-D 
method and the P-M method comparing to the original transfer method 
(Delta method). 

% 

In the sample problems reported in Appendix E, we have found that 
for straight, simple configurations (where the phenomenon of the columns 
becoming parallel is pronounced) the P-M and M-D methods give one or two 
more natural frequencies than can be obtained by the Delta method or 
the original transfer method. However, in more complicated configura- 
tions it seems to be possible to get three or four more natural fre- 
quencies. In the simpler cases, there are comparison values available 
from "exact" theory so that the accuracy of the results may be establish- 
ed. However, in the more complicated cases, such comparison values are 
not available. However, we have reasonable assurance of the integrity 
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of these results since we have been able to get the same values from 
the M-D method and from several different P-M subprocedures. Thus, 
using these methods, we seem to be able to double or triple the fre- 
quency range in which valid results may be obtained. 

Details of the systems analyzed and the results obtained may be 
found in Appendix E. 

3-2 Limitations and some suggestions. 

Non-stiff hangers may be introduced into the system by a point 
matrix with appreciable linear and/or torsional spring compliances. 

But, for stiff hangers the transfer method fails, so we cannot use 
program VIPIPE for natural frequencies of a system which has stiff 
hangers. In particular, the suggestion of Fink (4) that rigid hangers 
can be successfully handled by use of idealized very stiff exterior 
springs proves to be untenable. 

Marguerre and Uhrig suggested possible way to overcome the difficulty 
relating to stiff hangers by introducing unknowns [2] . Introducing un- 
knowns means method abandoning the transfer concept, so we cannot use 
VIPIPE directly to solve such a problem. However, the writer thinks 
that there is a possibility of modifying VIPIPE to solve such problems 
by getting state quantities at each node and assuming values for the 
unknowns. For details refer to reference (2] . 
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APPENDIX A 



DESCRIPTION OF PROGRAM 



A-l General remarks 

Program VIPIPE is a CDC FORTRAN 63 language digital computer program 
designed for use in the vibration analysis of planar piping systems, 
originally developed by George E. Fink (4j . The program VIPIPE is 
modified and augmented by the writer to use Delta method (same results 
as original transfer method, but with less amount of calculation), Modi- 
fied Delta method, and Pestel and Mahrenholtz' s remainder method with 
double precision or single precision arithmetic. 

In- and out-of-plane undamped natural frequencies may be found 
through the use of program VIPIPE. How this program may be used is 
fully explained in Appendix B. Considerable flexibility with regard to 
the mathematical model employed is provided in the program. For in- 
stance, one may use combination of a distributed mass treatment for 
straight section with a lumped mass treatment for curved sections (this 
method gives the greatest accuracy coupled with least machine time) or 
a lumped mass for all sections. Also one may consider or neglect the 
effect of shear deflection and rotational inertia, singly or together. 

Codes for the boundary conditions and various methods are inter- 
preted within the program. 

Problem solutions proceed by iterative processes with radian fre- 
quency the independent variable. The sign of the frequency determinant 
or remainder is used to control the search for and convergence to a 
system natural frequency. 
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Output and transfer to the next method/or transfer to the next pro- 
blem or to the end is provided in the program when it reaches signifi- 
cant figure limitation or previously set iteration limit, and also when 
a method fails. 

All subroutines are provided in binary deck in single precision 
and double precision arithmetic in the FORTRAN 63 version. Modifications 
of the program may be accomplished by changing FORTRAN cards 

only in the main program; the subroutines are in binary to conserve 
space and time. 

A-2 Program structure 

The program consists of a six section main body and twenty-seven 
subroutines. In the FORTRAN 63 library, GRAPH subroutine is not avail- 
able. So GRAPH subroutine is provided in binary, and built into the 
program. Also, in the FORTRAN 63 library, hyperbolic functions are not 
available, so these functions are computed in exponential form. 

The function of each section of the main body and of each subroutine 
is given below. 

A. Main body sections and their functions. 

INPUT Controls read-in of data, allocates storage loca- 

tions for arrays and subscripted variables, de- 
cides the first method to perform. 

INVARIANTS Computes and assigns appropriate subscripts to such 

system component characteristics as moment of in- 
ertia and radius of gyration, computes frequency 
increment for the first two modes, and sets the solu- 
tion acceptability criterion. 
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CONTROL 
IN -PLANE 


Constructs transfer matrix of the system for the in- 
plane case by calling for various subroutines in a 
sequence governed by the system geometry, the mathe- 
matical model specified, and the desired assumptions. 


CONTROL 
OUT-OF -PLANE 


Performs the same function as CONTROL IN-PLANE but 
for the out-of-plane case. 


ITERATION 


Evaluates the frequency determinant (modified fre- 
quency determinant in the M-D method) or the re- 
mainder and utilizes its sign in controlling the 
iterative search for mode frequencies. Initiates 
output when the specified frequency range has been 
fully explored, the required number of modes found, 
the specified iteration number has been reached, or 
the significant figure limit of computer reached. 


OUTPUT 


Controls output format, provides graph output, seeks 
the method to perform. 



B. Subroutines and their functions. 



SUBSEC 


Subdivides components to be treated by lumping mass 
on the basis of L/D ratio and starting frequency or 
required mode numbers. (L/D ratio for curved section 
means ratio of radius of curvature to diameter and 
included angle of arc, as applicable.) 


MATMUL 


Constructs the transfer matrix of those components 
treated by lumping mass for only one point mass and 
massless subsection. 


FINMAT 


Constructs the system state matrix (6x3) by successive 
premultiplication of the partial system state matrix 
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FINBRA 



STATEM 



STATEB 



BRCOR 



DELMA 



by the transfer matrix of each component as soon as 
each of the latter is constructed. For the lumped 
mass system, gives successive premultiplication by 
the transfer matrix constructed by MATMUL as many 
times as the number of lumped masses. 

Constructs the state matrix of a branch system by 
successive pre-multiplication of the partial branch 
state matrix by the transfer matrix of each component 
of a branch as soon as each of the latter is construct- 
ed. For the lumped mass system, gives successive pre- 
multiplication by the transfer matrix constructed by 
MATMUL as many times as the number of lumped masses. 
Constructs starting state matrix from the starting 
boundary condition of main member, for each method. 
Constructs starting state matrix of a branch system 
from the far end boundary condition of a branch for 
each method. Computes the purging factors of the 
branch in the M-D method. 

Constructs branch correction matrix for the assumed 
branch starting state matrix at a certain frequency 
in the P-M method. Constructs branch purging matrix 
in the M-D method. 

Constructs frequency determinant in the Delta method. 
Evaluates the purging factors from the frequency 
determinant of the Delta method and constructs the 
modified frequency determinant after another round 
of calculation in the M-D method. For the P-M 
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method, computes the correction factors for each sub- 
procedure. For each sub -procedure, if it fails by 
the fact that denominator is zero, the rows are 
changed and a second test is made. It also provides 
transferring to the next method or sub-procedure if 
it fails. 

DISTM, DISTMO, SFIELD, SFIELO, CFIELD, CFIELO, POINT, POINO, RIGID, RIGIO, 
STIFCO, HANGER, HANGEO, BRANCH, BRANCO, STAVEC , STAVEO, INVERT. 

For the above subroutines refer to reference [4] . 



A-3 



Program 

B1,B2,B3 



BC 



BL 

Dl,D2 



DEL 



DV 



EDL 



ERM 



FF 



nomenclature. 

Correction factors for assumed state quantities of 
the far end of a branch. 

Discriminant. Controls boundary correction in the P-M 
and H-D method. 

Number of the last iteration. 

Purging factors in the M-D method. Correction factors 
in the P-M method. 

Array of numerical values of frequency determinants 
and array of remainders.. 

Numerical value of a determinant in denominator of 
XI and X2. 

Array of numerical values of frequency determinants 
or remainders in single precision for input to GRAPH. 

A discriminant same as REM in the problem statement 
data card. 

Discriminant. Controls output to provide a print 



42 



FINK 



HH 



HMH 



K 

KA 



KKK 



MAT 



NUMPTS 
NUMEND 
PI, P2 
PMM 



output that analysis is terminanted due to signifi- 
cant figure limitation in the P-M method. 
Discriminant. Controls constructing the frequency 
determinant or modified frequency determinant in the 
M-D method. 

Discriminant. Controls transferring to another sub- 
procedure or to another way of testing in case one 
fails in the P-M method. 

Discriminant. Controls subsectioning in lumped mass 
system. 

Discriminant. Controls grid in the GRAPH output. 
Discriminant. Terminates iteration when specified 
iteration number is finished. 

Discriminant. Controls transferring to the output 
in case the last sub-procedure in the P-M method 
fails or controls transferring to another three sub- 
procedures (D,E,F) after three sub-procedures (A,B,C) 
are finished in case subprocedure C fails in Selec- 
tion 3 (for Selection, see Appendix B-9-6). 
Discriminant. Controls matrix multiplication in the 
lumped mass system. 

Number of points in one graph (except the last). 
Number of points in the last graph. 

Assumed starting state quantities of the main member. 
Discriminant. Controls three methods (Delta, M-D, 

P-M method) specified. 
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PPM Discriminant. Assures that same procedures are used 

for out-of-plane solution that were used for in-plane 
solution. 

PM Discriminant. Controls three methods to be used. 

PP Discriminant. Controls the three sub -procedures 

(A,B,C) in the P-M method. 

RR Discriminant. Controls testing by interchanging two 

rows in case one sub-procedure fails because the de- 
nominator of XI or/and X2 is zero in the P-M method. 

R3 Array name. Matrix product [rJ \ , used in 

the branch correction. 

REM Discriminant. Controls three ways of selecting sub- 

procedures in the P-M method. 

Ql,Q2,Q3 Assumed starting state quantities of a branch. 

SS Discriminant. Controls the three sub-procedures (D, 

E,F) in the P-M method. 

XX Array name. Branch correction matrix (3x3) in the 

P-M method and branch purging matrix in the M-D 
method. 

X1,X2,Y1,Y2 Correction factors for assumed starting state quanti- 
ties of the main member in the P-M method. 

YY Array name. Branch correction matrix (3x3x23). 

For other nomenclature see reference [ 4 ] 
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APPENDIX B 



INSTRUCTION FOR PROGRAM USE 

B-l General remarks. 

Appendix B is intended to provide complete instructions for using 
program VIPIPE. VIPIPE is FORTRAN 63 language program with twenty-seven 
subroutines. All subroutines are provided in a binary deck so that the 

compiling time is reduced as is the volume of the deck. 

< 

Program VIPIPE can be used in single precision or double precision 
arithmetic. Also VIPIPE is provided in FORTRAN 60 language with single 
precision. 

The Delta method, the M-D method, and the P-M method are all incor- 
porated in program VIPIPE. The M-D method and the P-M method do appear 
to permit obtaining more natural frequencies (one, two, or possibly several 
more) than can be obtained with the Delta method (same result as the 
original transfer method). However, the P-M method has the disadvantage 
of introducing an odd behavior of the remainder. In some cases, this 
odd behavior can easily be misinterpreted so as to lead to false fre- 
quency determination. However, there are a whole group of associated 
P-M modifications* and the odd behaviors occur at different frequencies 
for the different modifications. Thus, by making the calculation more 
than once, using a different modification each time, the true natural 
frequencies can be sorted out. 

Program VIPIPE uses the Delta method as the basic procedure. If 
more frequencies are required than have been obtained by this method, 

*P-M modifications means sub-procedures of the P-M method. 
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one may then use one or more of the available associated P-M modifica- 
tions and/or the M-D method. Graphical output will assist in identifying 
the natural frequencies. The writer recommends to use the P-M modifica- 
tions as original method in getting more natural frequencies, and to 
use the M-D method for purposes of comparison. 

The accuracy and number of natural frequencies that can be obtained 
using the methods depends on the capacity -of handling numbers accurately. 
As we can see in Appendix E, double precision arithmetic gives better 
results especially for complicated systems for which many matrix multi- 
plication are performed. But using double precision arithmetic results 
in considerable increase of machine time. And also for various methods, 
the machine time is different. The Delta method and the P-M method each 
take about the same amount of machine time, but the M-D method takes 
almost twice this time. 

So a user should consider the frequency range under investigation and 
the nature or characteristics of the system as well as the machine time. 

For comparison purposes Section B-9-6 of this Appendix lists machine 
times for one iteration for some example problems using each of three 
methods with double and single precision arithmetic. 

Complicated systems are divided into one main member and branches. 
System components are defined as sections of pipe, straight or of con- 
stant curvature, an elbow, a spring hanger, a coupling, a flange. Dia- 
meters, wall thickness, density, elastic and shear modulus may vary from 
component to component, but may not vary within a component. Both tors- 
ional and linear spring rates (for hangers) must have finite values. 

Thus perfectly rigid hangers cannot be treated. 
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B-2 Coordinate system. 

A local coordinate system is constructed at each point of interest 
in the system. A right handed coordinate system is used whose X axis 
is the local tangent to the central axis of the main member, always 
directed away from the left end. The X-Y plane coincides with the 
quiescent plane of system. The Z axis is always directed downward. The 
direction of the local Y axis is chosen so as to complete the right hand- 
ed system. 




Y 



Fig. B-2-1 

The branch coordinate starts from the end of the branch farthest from 
the joint point with the main member. The plane of the branch coincides 
with the X-Y plane. The positive Z axis of the branch and main member 
are parallel and have the same sense. The Y axis is oriented so as to 
complete a right handed system. 

B-3 Data required. 

Extensive properties. 

a. outside diameter (D) - inches. 

b. wall thickness (DI) - inches 

c. radius of curvature (RHO) - inches (Negative, if center of 
curvature lies on negative Y axis) . 
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d. length (TL) - inches (for straight section only; for 
curved sections, length is computed by the program). 

e. arc central angle (THETA) - degrees 

f. linear spring constant in X direction (CU() - pounds per 
inch 

g. linear spring constant in Y direction (CLY) - pounds per 
inch 

h. linear spring constant in Z direction (CLZ) - pounds per 
inch 

i. torsional spring constant about X axis (CTX) - in-lb/radian 

j. torsional spring constant about Y axis (CTY) - in-lb/radian 

k. torsional spring constant about Z axis (CTZ) - in-lb/radian 

Intensive properties 

3 

a. density of component (AAMU) - lb/ft 

b. elastic modulus (AE) - psi 

c. shear modulus (AG) - psi 

B-4 Component identification.* 

Components are identified by two numbers. 

A. All components have a component sequence number starting from 
the most left component of the main member. This sequence number is 
identified by the order of data cards containing component extensive and 
intensive properties. This can be understood from the figure B-4-1. 



*Aslo see pp. B-5 to B-7 of reference (4) . 
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B. A branch component identification number (NNBR) is read in 
as data. It indicates whether it is in a branch or main member, if it is 
in branch, the position in branch. 

NNBR 0 main member 

NNBR 1 the first component of a branch (farthest from the 
intersection) 

NNBR 2 intermediate component of a branch 
NNBR 3 the last component of a branch 

NNBR -1 the only component of a branch 

This can be understood from the figure B-4-2 




Fig. B-4-2 
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B-5 Branch joining angle 

Branch joining angle to main member is measured counter-clockwise 
in degrees from the X axis of the main member, as it is oriented immedi 
ately before the inter-section, to the X axis of the branch as it is 
oriented at the point of intersection. 



B-6 Boundary conditions* 

Five common boundary conditions, those of a fixed, free, pinned 
propped, or a roller supported end are indicated by the code given. 

The program interprets the code and forms the required state 
matrix (6x3). 



Boundary condition 


code: 


a. 


SV 1 


fixed end 


b. 


SV 2 


free end 


c. 


SV 3 


pinned end 


d. 


SV 4 


propped end 


e. 


SV 5 


roller supported end 


f. 


SV 0 


other than one of the above. 



(The user must construct the necessary state 
vector and cause to be read in as data) . 
This code is used for both in-plane and out-of-plane cases. 



> 



B-7 Control cards. 

A user sets up for execution with a master control, a FORTRAN con- 
trol card, and combination of END, FINIS, EXECUTE, transfer cards, bin- 
ary cards, and FORTRAN cards properly placed in FORTRAN 63 deck (JOB and 
END cards are used in FORTRAN 60). Figure B-7-1 shows the deck assembly 

*For further details, see pp. B-8 to B-10 of reference (4) . 
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for FORTRAN 63 use. 




B-8 Data deck assembly. 

The data deck is ordered as follows: 

a. problem statement card. 

b. boundary condition data card(s). 

c. branch joining angle data card(s). (deleted when system 
has no branches) 

d. component extensive properties data card(s). 

e. component intensive properties data card(s). 

The user may obtain solutions for any number of problems desired; 
in this case each problem data deck is assembled just behind the other. 
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Following the last problem deck a blank card should be added. 

B-9 Flexibility of program VIPIPE. 

B-9-1 Double precision and single precision arithmetic. 

Program VIPIPE can use double precision process only with FORTRAN 63. 
Double precision is used in system state matrix (UU) , branch system state 
matrix (W) , and frequency determinant or remainder. 

Type double is declared in main program and subroutines FINMAT, FIN- 
BRA, STATEM, STATEB, BRANCH, BRANCO, BRCOR, DEmA. By removing cards 
declaring type double the program can easily be changed into single pre- 
cision. Subroutines are provided in binary deck with single precision 
and double precision arithmetic, so with corresponding binary subroutine 
deck and main program in FORTRAN with or without a card (card number 
000120) declaring type double we can use single or double precision arith- 
metic. 

B-9-2 GRAPH. 

A user can get GRAPH output by option. The values of the frequency 
determinant and the values of the remainder are not usually of the same 
order of magnitude throughout the frequency range. Moreover, because 
of the infinite phenomenon of the P-M method, it is not general ly useful 
or convenient to get one graph for the whole frequency range. 

1. Number of points in one graph can be changed by changing 
NUMPTS in card 006090. 

2. Can get graph with grid or without grid by changing card 
006100. 

K = 1 with grid 

K = 0 without grid. 
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3. Graphs are labeled from 1,2, 3,... up to 10. Graphs beyond 
tenth one will oe labeled as 10. The last graph is label- 
ed LAST. 

4. For each graph, the whole frequency range is moved (except 
the first one) so that the first point of each graph starts 
from 0. The corresponding value of frequency determinant 
(or remainder) for the first point and NUMPTS (NUMEND for 
the last graph) are printed in the print output. 

B-9-3 Shear area form factor. 

A shear area form factor of one half is built in program.* This may 
be changed by changing card 000600. 

B-9-4 Starting frequency, frequency increment and solution acceptability 
criterion. 

Using the M-D or the P-M method, we are expecting to get more natural 
frequencies at high frequency. Thus choosing the proper starting fre- 
quency for the M-D or the P-M method, following the Delta method can 
save machine time. This can be done by putting the wanted starting fre- 
quency in problem statement data card. 

To compute the starting frequency increment, the program constructs 
a synthetic straight pipe, equal in length to the length of the main 
member, and having diameter, wall thickness, and intensive properties 
equal in magnitude to a weighted average of those properties of the 
composite members. One quarter of the fundamental frequency of this 
synthetic pipe is taken as the starting frequency increment (card 001390). 
And one thousandth of this increment is taken acceptability criterion 
(card 001410) . When the iteration reaches this solution acceptability 



*See reference (7) . 
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criterion, the corresponding natural frequency is determined using linear 
interpolation. 

This frequency increment is for the first two mode frequencies, but 
the acceptability criterion is the same for all mode frequencies. After 
second mode frequency is found, frequency increment is one tenth of the 
difference between the last two mode frequencies found (card 003950) . 

These can be altered by changing cards. The writer found that a 
smaller frequency increment gives better results especially at high fre- 
quency with the P-M method. 

B-9-5 Subsectioning. 

Subsectioning for the lumped mass treatment of components is carried 
on in subroutine SUBSGC. A component is divided into maximum of twelve 
subsections. 

For all components having a length to diameter ratio greater than 
six, the number of subsections is computed as two times the value of a 
parameter HM which appears in SUBSEC (HM is number of modes to be sought) 
If HM is greater than six, HM is given the value of six. In case the 
starting frequency is greater than 800 radians per second, HM is also 
given the value of six. HM can be controlled by a parameter HMH in the 
main program (card 001250). If HMH is zero, HMH does not affect the 
value of HM. By changing HMH to other than zero we can change the value 
of HM equal to that of HMH. 

Poor subsectioning might give some loss of accuracy in determining 
natural frequencies, especially at high frequency. Checking the natural 
frequencies by increasing the number of subsections using parameter HMH 
is recommended. 
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B-9-6 Machine time and printing line limit. 

In FORTRAN 63, estimated time is set to thirty minutes and print- 
ing line limit is set at 30000 lines in VIPIPE. These limits can be 
altered by changing the master control card. If a limit is reached, 
computation stops automatically. 



For reference, machine time for one iteration of some example problems 
is given below. 





Example System* 


Delta meth. 


M-D meth. 


P-M meth. 




system 


1 


0.131 sec. 


0.247 sec. 


0.144 sec. 


4) • 

r-» U 
Xi 4) 
n Li 


system 


4 


2.10 sec. 


4.22 sec. 


• 2.14 sec. 


O PM 
Q 


system 


5 


4.799 sec. 


6.90 sec. 


3.505 sec. 


4) • 

T-t O 

hA It 


system 


1 


0.044 sec. 


0.117 sec. 


0.045 sec. 


W V 

c u 
pm 
CO 


system 5 


0.377 sec. 


0.968 sec. 


0.49 sec. 



B-9-7 Performing various methods. 

One selects the methods to be employed by use of the problem state- 
ment data card. Then the program picks the first method in the INPUT 
section of the main program with a discriminant PM, in the order: Delta 
method, M-D method, P-M method. 

Finishing a method, the OUTPUT section of the main program opens to 
the next method. Finishing all specified methods, the program goes to the 
out-of-plane portion, to the next problem, or to the end. 

In the P-M method, there are six sub-procedures. These sub-pro- 
cedures are controlled by discriminatns PP, SS. 



PP : 


0. 


subprocedure A, 


SS : 0. 


sub-procedure D, 


PP : 


1. 


sub-procedure B, 


SS : 1. 


sub-procedure E, 


PP : 


2. 


sub-procedure C, 


SS : 2. 


sub-procedure F. 



*For example systems refer to Appendix E. 
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Each sub-procedure has another possible way by interchanging rows of final 
determinant (refer to section 2-4). If one fails, the alternative one 
is tested. If one succeeds, then go to the next sub -procedure. There 
are three ways of selecting sub -procedures by discriminant REM. 

REM: 2. three sub-procedures in the order: A,B,C. - Selection 1 

REM: 1. three sub-procedures in the order: D,E,F. - Selection 2 

REM: 0. six sub -procedures in the order: A,B,C,D,E,F. - Selection 3 

If one wants just one sub -procedure, by changing PP to a value of 2 
(card 000710) in selection 1 only sub-procedure C is used, or changing 
SS to a value of 2 (card 001000) in selection 2 only sub-procedure F is 
used. If one wants just two sub -procedures , by changing PP or SS to a 
value of 1 sub-procedure B and C or sub-procedure E and F, respectively, 
are used. 

B-9-8 Limiting the number of iterations. 

One can limit the number of iterations, for each problem, by a para- 
meter HL in the problem statement data card. In the P-M method, because 
of the infinite phenomenon, the frequency increment in iteration can be- 
come very small. Thus the iteration might keep going on indefinitely. 

The program permits six hundred iterations as maximum. The writer 
found three hundred Iterations to be adequate to get nine natural fre- 
quencies with the acceptability criterion set in the program. 

B-9-9 End conditions. 

A boundary condition state vector must have three non-zero and three 
zero elements. Frequently, however, an actual condition may not meet 
this specification. For instance, an end may be supported or connected 
to a piece of machinery. Then we approximate it as being held in a 
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mounting block having specified linear and torsional spring rates, and 
to this attach a fictious massless short section. 

Through the use of this sort of fictious end, a proper end condi- 
tion state vector may be constructed for any real situation with little 
or no loss of accuracy. 

For hangers exhibiting coupling, systems having unusual flexibility, 
and flanges and couplings refer to pp. B-19 to B-20 of reference (4) . 

B-10 Format of data 
B-10-1 Problem statement 



The card has 15 fields under control of the following field specifi- 
cations. 



312, 


4F3.0, 


2F8.0, 4F2.0, 2F4.0 


Field names and the information each field conveys are given below in 
order. 


A. Field 


names 


and content. 


1. 


NS 


number of components. 


2. 


NBR 


number of branches. 


3. 


IOP 


code for in-and/or out-of -plane solutions. 


4. 


HM 


number of modes sought. 


5. 


AMYK 


code calling for the use of lumped or distri- 
buted mass model 


6. 


SD 


code calling for shear deflection inclusion or 
neglect 


7. 


RI 


code calling for rotational inertia effect in- 
clusion or neglect 


8. 


OMEGAG 


upper limit of frequency range of investigation 


9. 


CMEGAI 


lower limit of frequency range of investigation 
and initial iteration frequency. If 0 is placed 
it will begin at .01 radians per sec. 



57 



10. 


GDABC 


code calling for print output. . 


11. 


BRC 


code calling for the inclusion or exclusion of 
branch correction factors in P-M method and 
purging factors in M-D method. 


12. 


GRAPH 


code specifying whether there is to be graph 
output. 


13. 


REM 


code for choosing methods in P-M method. 


14. 


PMM 


code calling for the combination of 3 methods of 
Delta, M-D, P-M method which a user wants. 


15. 


HL 


number of iterations. 


Codes 


used. 




1. 


IOP 0 


Both in-and out-of-plane mode frequency sought 
(in this case, in-plane frequencies are sought 
first) . 


2. 


IOP 1 


In-plane frequencies are sought. 


3. 


IOP 2 


Out-of -plane frequencies sought. 


4. 


AMYK 1. 


Lumped mass model. 


5. 


AMYK 0. 


Distributed mass model, 


6. 


SD 0. 


Consider shear deflection. 


7. 


SD 1. 


Neglect shear deflection. 


8. 


RI 0. 


Consider rotational inertia. 


9. 


RI 1. 


Neglect rotational inertia. 


10. 


GDABC 0. 


Print solution only. 


11. 


GDABC 1. 


Print solution and input data. 


12. 


GDABC 2. 


Print solution and input data and graph data. 


13. 


BRC 0. 


Do not consider branch correction factors or 
purging factors. 


14. 


BRC 1. 


Consider the branch correction factors or purg- 
ing factors. 


15. 


GRAPH 0. 


No graph output. 


16. 


GRAPH 1. 


Graph output. 



58 



17. 


REM 


0. 


All six sub-procedures 
C,D,E,F). 


in the P-M method (A,B, 


18. 


REM 


1. 


Three sub -procedures in 


the 


P-M method 


(D.E.F) 


19. 


REM 


2. 


Three sub-procedures in 


the 


P-M method 


(A,B,C) 


20. 


PMM 


1. 


Delta method only. 








21. 


PMM 


10. 


M-D method only. 








22. 


PMM 


100. 


P-M method only. 








23. 


PMM 


11. 


Delta, M-D method. 








• 

eg 


PMM 


101. 


Delta, P-M method. 








25. 


PMM 


110. 


M-D, P-M method. 








26. 


PMM 


111. 


Delta, M-D, P-M method. 









B-10-2 Boundary conditions. 

Field specification 22F7.3 

Boundary conditions are read in separately for in- and out-of-plane 
cases. When both cases are to be investigated, in-plane case data are 
read in first. 

N + 2 fields must be filled in (N is number of branches) . The first 
and last fields are first and last boundary condition codes of the main 
member. Branch boundary conditions are entered in the intermediate fields 
in accordance with the branch sequence number. When a boundary condition 
is not describable by one of the codes, a zero must be entered in its 
data field. For each zero an additional data card is required. The 
limitation in non-zero state quantitites at the boundary must be three 
(refer to B-9-9) . 

B-10-3 Branch joining angle (degrees). 

Field specification is: 20F7.3 
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B-10-4 Component extensive properties and branch component identifica- 
tion. 

A data card for each system component, ordered by sequence number. 
Field specification is: 3F5.2, 2F6.2, 6F7.2, 12 

Field names are: in order, 

D, DI, RHO, TL, THETA, CLX, CLY, CLZ, CTX, CTY, CTZ, NNBR 
B-10-5 Intensive properties. 

The first intensive properties data card has four fields with speci- 
fication: 

F7.3, 2E10.2, F2.0 

Field names are: AAMU, AE, AG, DD* 



*When the intensive properties of all components are the same, the DD 
field is left blank. When the intensive properties are different, the 
DD field is non-zerb, then an intensive data card must be prepared 
for each of the remaining components (including hangers) . Additional 
data cards are identical to the first one except that the DD field is 
blank. 
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APPENDIX C 
FLOW DIAGRAMS 



C-l General remarks. 

Flow diagrams are included for all parts of the main program, and 
for several subroutines. Since CONTROL OUT OF PLANE, FINBRA flow dia- 
grams are identical in structure to those of CONTROL IN PLANE and FINMAT 
respectively, only the latter are included. Reader can refer to the 
flow diagram of STAVEC and STAVEO of reference [4]. Since the flow 
diagram structures of SUBSEC, BRANCH and BRANCO are similar to those 
of reference [4], they are not included here. 

Flow diagrams of other subroutines that are not included in this 
paper are relatively simple ones, so the reader can refer directly to 
the program. 

C-2 Index of flow diagram. 



Main program. 

Flow Diagram. Page. 

INPUT 63 

INVARIANTS 65 

CONTROL IN PLANE 66 

ITERATION 68 

OUTPUT 70 

Subroutines 

FINMAT 73 

STATEM 74 

STATEB 75 

BRCOR 76 

DELMA 77 
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C-3 Flow Diagrams 



INPUT PLOW DIAGRAM 
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INPUT FLOW DIAGRAM(CONT. ) 





I •■.•VARIANTS FLOW DIAGRAM 



FROM I!?PUT 






/CLX(l)-CLY(l) = CLZ(l)\ 

‘ CTX(l)*CTY(l) = ctz(i) 



1= 1,NS 






AMP(l)aaAMU(l)«ny( D(l) »»2- Dl(l)»»2)/( 691 2. *32.17*12.) 



THETA(l)-= THETA(l)**/l80. 

tl(i)~theta(i)*absf(rho(i)) 


*%o U> 

JL =o 








) 


AJT( I )- 7t*(D( l)**4Dl( I ) **4)/32. 
AJ(l)=AJT(l)/2. 

r(i)= i./ ( aj(i)*e(i)) 

AIX(I)« SQRTF((B(l)**2 +Dl(l)**2)/8. 
AIY(I) =» AIX(l)/l. 414214 




> 


t 



MTS 




0 ( 1 ) 

=0 









^ — > 


HHM=RM 0MH=O. 


■ CALL 
KM= 


SUBSEC 

HHM 









AL=-AL + TL(l) 

AD =■ AD+ B(I)*TL(I) 
EE=rEE-|-E(l)*TL(l) 
$[AID=sAID+Dl(l)*TL(l) 
AAMU— AAMU f- AMU( I ) *TL( I ) 




* 


£0 

t .... 


AD - AD/AL 
AID = AID/AL 
AFR »T?64.*(AD* 


, 1 
*4-AID**4) 


ADOM *» FR/4.0 
DOM - ADOM 
CR= FR*0. 00025 




FR = 4.73**2»SQR?F(EE*AFR/(AAMU*AL**4) 


BDOM - ADOM 


BR= CR 





5 
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INVARIANTS r L OY P I AGRAM( CO NT. ) 




CONTROL IN PLANE FLOW DIAGRAM 




CONTROL IN PLANE DIAGRAM (CONT.) 





ITERATION FLOW DIAGRAM 







ITSRATIOH FLOW DIAGRAM(CONT . ) 



delc«del(l) 

OMEGAC =OMEGA(L) e 

~JT" 




— <^A>f (bB^oV> — ^<PDEL*DEL( ] 

! 0 

40/ ^ >0 



<Q 



(a)^— |'omega( l~m ) = omegaClKdom k— — < oh¥ga( l Pom e g7g> ^ <aa>— — 3 



I 0MEGA ( Lfl ) =■ om ega(l ] 



-(omega(l)-omegac)/2. 



A 



(£V - |cc=rl7V ^ /0MEGA(L)-0M EGA(L-l )-CRy 




OMEGAC ** OMEGA(L-l) 
AA = 1. • 

BB=1. 

DELC=i)EL(L-l) 



o k bg am(m)= o„ E o,c, 



BM*mM 






>o 



>o 




:-BM>H K= Mtl 

±0 




XL 



XO 



ADOM«= ( 0MEGAM(M-l)-0MEGAM(M-2) )*. 1 



DELC = 0 . 

AA~BB-CC =GG=0 

DOM - ADOM 

7E 



It- 



TO OPT PUT 



IoMEGA( L-*~l)— 0MEGAM(M-l)f‘2».*CR 



GO i58o Address 

CONTROL IN PLANE 



OMEGAO IOMEGA ( L~tl )[~>(A) 



(k) — j |L » LtlW iL 




iTO OUTPUT^ \ KA=» 1 | 







OMEGA(lrfl) • 0MEGA(L)-.9*D0M 
DOMsr . 5*DOM 



GO TO Addressl002 
CONTROL OUT OP PLANE 



6 ? 



OUTPUT FLOW DIAGRAM 




TO 









OU TPUT FLQV: PI AGRAM(CONT» ) 
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OUTPU T PLOW DI AGR AM( CONT . X 




'll 




SUBROUT INE FINM AT FLOW DIAGRAM . 
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SUBR O UTINE STAT EM FL OW DIAC RAM 




74 - 





SUBROUTINE ST AT OB FLOW DIAGRAM 




|v( J,K)“YY( j,k,n) 1 
1 J = l,3 K^1.3 I 



CALCULATE: 
Purging factors 
D1 & D2 



, W(J,K)=0. 

1 



W(J,2) = -D1 
VV( J, 3) =, -D2 




Bl— B2-B3-0. 
Dl= D2=0. 



= 0 



<0 






^0 


hTn 

> 

| 


) — YY( J ,k,n) 
,3 K-1,3 








-5T- 

CALCULATE: 
B1,B2,B3. 

For subproc. 
A &D 



CALCULATE: 
B1,B2,B3. 
For subproc. 
C & F 



CALCULATE: 
B1,B2,B3. 
For subproc. 
B & E 



I 



Q1~Q1*B1 
Q2=Q2*B1 B2^ 
Q3= *Q3* B1 B5 



—0 



|M=1 




jf.. 






& 
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SUBRO UTI NE STAV EB ■ .LOW. DIAGRAM( CONT. 






W(J,1)-Q3 

vv(j,m)-i. 



>0 



■±. 



>0 




& 



<0 



T =*& 


i 


>0 

t 3 


i 


W(J,1)-Q2 


vv(j ,l)=Ql 




w(j,m)=i. 










1 , 




\/ 






M-Hl-ll 
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APPENDIX D 



PROGRAM LISTING 



Sequence Number 



Page 



D-l INDEX 

Main Program. 

Section 
INPUT 

INVARIANTS 
CONTROL IN PLANE 
CONTROL OUT OF PLANE 
ITERATION 
OUTPUT 

i. \ 

Subroutines 

V iV 

SUBSEC 

DISTM 

DISTMO 

RIGID 

RIGIO 

STIFCO 

STIFOO 

STAVEC 

STAVEO 

INVERT 

HANGER 

HANGEO 

CFIELO 

POINO 



000000 - 001070 82 

001080 - 001500 83 

001510 - 002450 84 

002460 - 003390 86 

003400 - 004420 88 

004430 - 007990 89 

008000 - 008520 96 

008530 - 009090 97 

009100 - 009660 98 

009670 - 009920 99 

009930 - 010180 99 

010190 - 010430 100 

010440 - 010680 100 

010690 - 010990 101 

011000 - 011260 101 

011270 - 012000 102 

012010 - 012180 103 

012190 - 012360 103 

012370 - 012880 104 

012890 - 013100 105 
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Section 


Sequence Number 


Page 


CFIELD 


013110 


- 


013560 


105 


POINT 


013570 


- 


013790 


106 


SFIELD 


013800 




014150 


106 


SFIELO 


014160 


- 


014420 


107 


MATMUL 


014430 


- 


014560 


107 


FINBRA 


014570 


- 


014850 


108 


FINMAT 


014860 


- 


015140 


108 


BRANCH 


015150 


- 


015780 


109 


BRANCO 


015790 


- 


016470 


110 


STATEM 


016480 


- 


017020 


111 


STATEB 


017030 


- 


017780 


112 


BRCOR 


017790 


- 


018100 


113 


DELMA 


018110 




019500 


114 



D-2 Listing. 
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PROGRAM VIPIPE 



INPUT 



DIMENSION AJT( 50) ,AJ< 50) .R (50) . AIX( 50) .AI Y< 50 ) .PHI < 20) .SVI (6.1 ) » 
1DT( 50) » D( 50 ) » D I ( 50 ) , TL < 50 ) .Z ( 5" ) .G ( 50 ( , E ( 50 > . AMUI 50 ) .THETA ( 50) ,RHO 
2(50) .OMEGAMI 20) .DEL ( 600) .OMEGA ( 600 ).A(6.6).B(6»6)»U(6.6).CLX(50).C 
3LY( 50) . CLZ ( 50) .CTX(50) .CTY< 50 ) , CTZ ( 50 ) , V ( 6 . 3 ) . VV ( 6 . 3 ) »UU( 6 . 3 ) .SVE< 
46. 1 ) »SV( 22 ) .SVI P(6 .22 ) .SV0( 22 ) .SV0P( 6.22 ) .NNBR ( 50 ) .BAMU ( 50 ) . THETR ( 
520) .P( 20) »YY< 3.3.23) »R3 < 3.3) .XX( 3.3) » Inti-El 12) .E0L( 600) 

TYPE DOUBLE VV .UU . V. X2 . XI . Y2 » Y 1 .PI »P2 . D1 .D2 . DEL 
NM = 0 
JJ=1 ' 

450 READ 1500.NS.NBR. IOP.HM.AMYIC.SD.RI .OMEGAG .OMEGA I , GDABC . BRC .GRAPH. 
1REM.PMM.HL 

1500 FORMAT) 312. 4F3. 0. 2F8 .0 .4F2 .0 » 2F4. 0 ) 

IF (NS) 600.600.601 

601 NMBR=NBR+2 

IF( IOP-D 1501. 1501. 1506 

1501 READ 1502»(SV( I ).I=1.NMBR) 

1502 FORMAT ( 22F3.0 ) 

DO 1505 J= 1 . NMBR 

I F ( SV ( J ) 11503.1503.1505 

1503 READ 1504, ( SVIPI I . J) . I »1.6 ) 

1504 FORMAT ( 6F2 . 0 ) 

1505 CONTINUE 

IF(IOP) 1506,1506.1523 

1506 READ 1 507 , ( SVO (I), 1=1. NMBR ) 

1507 FORMAT ( 22F3 .0 ) 

DO 1510 J* 1 ♦ NMBR 

I F ( SVO( J) 11508,1508,1510 

1508 READ 1509»(SVOP(I»J)»I=1,6( 

1509 FORMAT ( 6F2 • 0 ) 

1510 CONTINUE 

1523 IF(NBR) 1513.1513,1511 

1511 READ 1512, (PHI(I). 1 = 1, NBR) 

1512 FORMAT (10F7.3) 

1 5 1 30READ 1514, (D(I), DKI ) ,RHO( I ) » TL ( I ) .THETA ( I ) ,ClX( I ) ,CLY< I ) ,CLZ( I ) »C 
1TX( I ) ,CTY( I ) ,CTZ( I) .NNBR ( I ) ,1 = 1. NS) 

1514 FORMAT! 3F5.2.2F6.2 .6F7.2.I2 ) 

READ 1515 , AAMU.AE.AG.DD 

1515 FORMAT ( F7 . 3 . 2E 10. 2 , F2 . 0 ) 

IF(DD) 1516,1516.1518 

1516 DO 1517 1=1, NS 
AMU (I )=AAMU 

E ( I ) = AE 

1517 G ( I ) = AG 

GO TO 1520 

1518 AMU ( 1 ) = AAMU 
E( 1 )=AE 

G ( 1 ) =AG 

READ 1519, (AMUI I ) , E ( I ) ,G ( I ) , I =2 .NS ) 

1519 FORMAT ( F7. 3 , 2E 10 . 2 ) 

1520 I F ( OMEGA I ) 1121.1121,1521 



000000 

000010 

000020 

000030 

000040 

000050 

000060 

000070 

000080 

000090 

000100 

000110 

000120 

000130 

000140 

000150 

000160 

000170 

000180 

000190 

000200 

000210 

000220 

000230 

000240 

000250 

000260 

000270 

000280 

000290 

000300 

000310 

000320 

000330 

000340 

000350 

000360 

000370 

000380 

000390 

000400 

000410 

000420 

000430 

000440 

000450 

000460 

000470 

000480 

000490 

000500 

000510 

000520 

000530 

000540 

000550 



82 



U \J U KJ 



1121 OMEGA I = *0 1 

1521 OMEGA( 1 ).®OMFGAI 
AMYK = AMYK.-1 . 0 
OME GAO® OMEGA I 
FFAC**5 

L* 1 

M s 1 

MN*0 

AA=0.0 

BB®0 *0 

CC=0.0 

FF*0.0 

GG = 0 • 0 

RRrO.O 

KK*0 

PP*0. 

KA = 0 
KKK*0 X 
AAMU*0 • 0 
EE*0 *0 
AL x 0 • 0 
AD®0*0 
NNN*NBR+1 

I F ( PMM*0# 1—1 • ) 6091*6092*6093 

6091 PM=—1 • 

GO TO 7000 

6092 PM=0 • 

GO TO 7000 

6093 I F ( PMM*0# 0 1-1 • ) 6094*6095*6096 

6094 PM*-1* 

GO TO 7000 

6095 PM* 1 • 

GO TO 7000 

6096 IF(PMM#0. 01-1*1) 6097*6098*6099 

6097 PM=-1 • 

GO TO 7000 

6098 PM® 0# 

GO TO 7000 

6099 PM*- 1 • 

7000 PPM*PM 

IF ( PM) 7073*7073*7001 

7001 I F ( REM— 1 • ) 7071*7072*7071 

7071 SS*-1. 

GO TO 7073 

7072 SS=0. 

7073 DO 1522 I*1*NS 
DT ( I ) *D I ( I ) 

BAMU ( I ) * AMU ( I ) 

THETR ( I ) *THETA( I ) 

1522 Dt ( I ) =D ( I )-2*0*DI ( I ) 



INVARIANTS 
DO 108 I * 1 * NS 

I F ( CLX ( I ) +CLY ( I )+CL2( I )+CTX( I )+CTY< I) +CTZ ( I ) ) 108 * 10 1 • 10 8 



000560 

000570 

000580 

000590 

000600 

000610 

000620 

000630 

000640 

000650 

000660 

000670 

000680 

000690 

000700 

000710 

000720 

000730 

000740 

000750 

000760 

000770 

000780 

000790 

000800 

000810 

000820 

000830 

000840 

000850 

000860 

000870 

000880 

000890 

000900 

000910 

000920 

000930 

000940 

000950 

000960 

000970 

000980 

000990 

001000 

001010 

001020 

001030 

001040 

001050 

001060 

001070 

001080 

001090 

001100 

001110 



1 

-t. 



S3 



nnnn 



10 1 AMU ( 1 ) =AMU ( I )*3.141592654*<D< I ) **2-DI (I)**2)/<6912.*32.17*12.) 
I F ( RHO ( I ) ) 102.103.102 

102 THETA ( I ) =THETA ( I ) *3 • 14 1 592654/ 1 80 • 

TL( I ) = THE T A < I )*ABSF< RH0( I ) ) 

103 AJT ( I )=3.141592654*(D( I )**4-Dl ( I )**4)/32. 

AJ ( I ) = AJT ( I ) / 2 • 

R ( I ) = 1 .0/ ( AJ < I ) *E ( I ) ) 

A I X .( I >=SQRTF( (D( I ) **2+DI ( I ) **2 ) /8. > 

A I Y ( I ) *A I X ( I )/l. 414214 
IF(AMYK) 104.105.104 

104 Z ( I )*1.0 

I F ( RHO ( I ) ) 105.106.105 

105 HHM a HM 
HMH=0. 

CALL SUBSEC ( RHO< I ) . THETA ( I ) . D< I ) . TL ( I > .HM .Z (I) .OMEGA I .HMH ) 
HM*HHM 

106 1F(NNBR( I )) 108. 107.108 

107 AL=AL+TL ( ! ) 

AD=AD+D ( ! )*TL( I ) 

EE = EE+E ( I ) *TL C I ) 

A I D = A I D+D I ( I ) *TL ( I ) 

AAMU C AAMU+AMU ( I )*TL( I ) 

108 CONTINUE 
AD=AD/AL 

A I D = A I D/AL 

AFR*3. 14159 2654/64.* ( AD**4-A 1 0**4 ) 
FR=4.73004019**2*S0RTF(EE»AFR/( AAMU*AL**4) ) 

AD0M=FR/4 .0 
DOM* ADOM 
CR = AD0M*0 *001 
BDOM=ADOM 
BR = CR 
ERM*REM 

DO 109 I ! * 1 • NBR 

109 PH HI I ) =PH I ( I I ) *3 *14 15 926 54/ 180.0 
I F ( IOP-1 ) 110.999.998 

110 MN r 1 



CONTROL IN PLANE 

999 P1M.0 
P2* 1 .0 
DUO. 

D2 = 0. 

Ql* 1 .0 
02 = 1.0 
03=1.0 
B1 = 0. 

B2 = 0. 

B3 = 0 • 

F I NK = 0 . 

IF(SV( 1) ) 1011.1021.1011 
1011 CALL STAVEC ( SV ( 1 ) . 1 . S V I ) 
GO TO 1041 
1021 DO 1031 J = 1 .6 



001120 
001130 
001140 
001150 
001160 
001170 
001180 
001 190 
001200 
001210 
001220 
001230 
001240 
001250 
001260 
001270 
001280 
001290 
001300 
001310 
001320 
001330 
001340 
001350 
001360 
001370 
001380 
001390 
001400 
001410 
001420 
001430 
001440 
001450 
001460 
001470 
001480 
001490 
001500 
001510 
001520 
001530 
00154Q 
001550 
001560 
001570 
001580 
001590 
001600 
001610 
001620 
001630 
001640 
001650 
001660 
001670 



S4- 



1031 SVI < J * 1 )*SVIP<J.l) 

1041 IF(SV(NMBR) >1051*1061*1051 
1051 CALL STAVEC ( SV ( NMBR ) t 1 *SVE ) 

GO TO 1081 
1061 DO 1071 J = 1 * 6 
1071 SVE ( J * 1 ) =SV I P ( J *NMBR ) 

1001 JK=0 

IF(NBR) 1000 * 1000* 1085 

1085 DO 1087 N=2*NNN • 

IFCSV(N) ) 1087*1087*1086 

1086 CALL ST AVEC ( SV ( N ) . N* SV I P ) 

1087 CONTINUE 
1000 N-2 

11*1 
HH*0 • 

DO 1001 I=1*NS 
MAT = 0 
BC*0 • 

P( I )=CLX< I ) +CLY ( I ) +CTZ ( I ) 

IF(P( I) >37*38*37 

37 CALL HANGER (CLX( I ) *CLY( I ) *CTZ ( I > *U) 

GO TO 62 

38 I F ( AMY 0 39*41 *39 

39 I F ( RHO ( I ) >40*42*40 

40 IF(Z( I) ) 43*45*43 

41 I F ( RHO ( I ) >40*48*40 

42 I F C Z ( I ) >46*47*46 

43 CALL CF I ELD ( R ( I ) *Z < I ) *G( I ) *SD*RHO< I ) ♦ THETA< I ) *D< I ) * DI ( I> • E < I ) * A ) 

44 CALL POINT ( AMU ( I ) *A I Y ( I ) *OMEGA( L ) *TL ( I > *Z( I > *R I *B) 

CALL MATMUL <A*B*U) 

MAT = 1 
GO TO 62 

45 CALL STIFCO ( THETA ( I > *RHO ( I ) *U ) 

GO TO 62 

460CALL DISTM <AMU( I > * A I Y ( I > * OMEGA ( L > * T L ( I > *RU ) *D< I > *DI ( I > *G( I ) *SD*R 
II *E( I > * FFAC *U ) 

GO TO 62 

47 CALL RIGID ( AMU ( I ) *TL ( I ) * OMEGA (L>*AIY(I)*RI*U) 

GO TO 62 

48 IF ( Z ( I > >49*47*49 

49 CALL SF I ELD ( DI ( I ) *TL ( I ) *R ( I ) *Z I I > *G( I > *SD*D( I ) *E ( I > *FFAC* A ) 

GO TO 44 

62 IF ( NNBR ( I ) ) 81*80*81 
80 I F ( 1-1 >63*79*63 

79 CALL STATFM(SVI*UU*D1*D2*P1*P2*PP*PM*SS) 

GO TO 63 

63 CALL F I NMAT (U*UU*V*MAT *A*Z( I > > 

DO 64 J*1 *6 

DO 64 K* 1 * 3 

64 UU< J*K)=V( J*K> 

IF(PM) 88*1113*84 

1113 I F ( FI NIC) 88*84*08 
84 IF(BC) 88*88*1111 
1111 I F ( BRC > 88*88*89 

89 CALL BRCOR ( R 3 *UU * XX * JK * VV * PM ) 

MMM=N-1 



001680 

001690 

001700 

001710 

001720 

001730 

001740 

001750 

001760 

001770 

001780 

001790 

001800 

001810 

001820 

001830 

001840 

001850 

001860 

001870 

001880 

001890 

001900 

001910 

001920 

001930 

001940 

001950 

001960 

001970 

001980 

001990 

002000 

002010 

002020 

002030 

002040 

002050 

002060 

002070 

002080 

002090 

002100 

002110 

002120 

002130 

002140 

002150 

002160 

002170 

002180 

002190 

002200 

002210 

002220 

002230 



25 



n on n 



DO 90 J=1 .3 
DO 90 K=1 .3 

90 YY( J.K.MMM) =XX( J.K) 

88 GO TO 1001 

81 IF ( NNBR ( I )-l ) 83.83.65 

83 CALL STATEBlSVIP.VV.N.YY.L.BRC.Dl. 

65 CALL FINBRAIU.VV.V.MAT .A.ZII) ) 

DO 66 J = 1 .6 

DO 66 K=1 .3 

66 VV< J.K)=V< J.K) 

I F ( NNBR ( I I 168.67.67 

67 IF(NNBR(I)-3) 1001,68.68 

68 CALL BRANCH(R3.W»PHI ( I I I »U) 

N=N+1 

11 = 11+1 
MAT = 0 
BC* 1 • 

GO TO 63 
1001 CONTINUE 

GO TO 2000 



CONTROL OUT OF PLANE 

998 Pl=1.0 
P2 = 1 .0 
D 1 = 0 • 

D2 = 0. 

Ql = 1 • 0 
02=1.0 
Q3= 1 .0 
81 = 0. 

B2 = 0. 

83 = 0. 

F I NK=0. 

IF(SVOIl) 11013,1023,1013 
1013 CALL ST AVEO ( SVO ( 1 ) .l.SVI ) 

GO TO 1043 
1023 DO 1033 J= 1 .6 
1033 SVHJ.l)=SVOP<J.l) 

1043 IFISVO(NMBR) I 1053.1063.1053 
1053 CALL STAVEO(SVOINMBR) .l.SVE) 
GO TO 1083 
1063 DO 1073 J= 1 , 6 
1073 SVE ( J» 1 )=SVOP (J.NMBR) 

1083 JK= 1 

IF(NBR) 1002,1002,1093 

1093 DO 1095 N*2 »NNN 

I F ( SVO ( N I ) 109 5.1095,1094 

1094 CALL STAVEO ( SVO ( N ) .N.SVOP) 

1095 CONTINUE 
1002 N=2 

11=1 
HH = 0. 

DO 1003 1=1, NS 
BC = 0. 



002240 

002250 

002260 

002270 

002280 

,PP,PM,SS.Q1»Q2.Q3.FINK) 002290 

002300 

002310 

002320 

002330 

002340 

002350 

002360 

002370 

002380 

002390 

002400 

002410 

002420 

002430 

002440 

002450 

002460 

002470 

002480 

002490 

002500 

002510 

002520 

002530 

002540 

002550 

002560 

002570 

002580 

002590 

002600 

002610 

002620 

002630 

002640 

002650 

002660 

002670 

002680 

002690 

002700 

002710 

002720 

002730 

002740 

002750 

002760 

002770 

002780 

002790 



86 



MAT = 0 

P( I ) =CTX ( I ) +CLZ ( I )+CTY( I ) 

IF(P( I ) ) 1*2*1 

1 CALL HANGEO ( CTX( I ) *CLZ ( I ) *CTY( I)*U) 

GO TO 14 

2 I F ( AMYK ) 3 • 1 1 * 3 

3 I F ( RKO ( I ) ) 4 * 8 *4 

4 I F ( Z ( I ) ) 6 . 5 * 6 

5 CALL STIFOO < THET A < I ) , RHO < I ) ,U ) 

GO TO 14 

60CALL CF I ELO ( A JT ( I).R(I)*Z(I)*G(I) *SD.RHO( I ) • THETA ( I ) * D ( I ) • D I ( I ) * FF 
1 AC • A ) 

7 CALL POINO ( AMU (I)*AIX(I)*AIY(I) * OMEGA ( L ) • TL ( I ) *Z ( I ) * R I * B ) ■ 

CALL MATMUL <A*B*U) 

MAT = 1 
GO TO 14 

8 I F ( Z ( I ) ) 10.9*10 

9 CALL RIGIO ( AMU ( I) * T L ( I ) * A I X ( I ) .OMEGA ( L ) . A I Y ( I) . A JT (I ) * R I . U ) 

GO TO 14 

100CALL D I STMO ( AMU ( I).AIX(I)*AIY(I) *OMEGA ( L ) *TL ( I ) * A JT ( I ) . R < I ) , D ( I ) . D 
1 1 ( I ) *G( I ) .SD.RI . FFAC *U ) 

GO TO 14 

11 I F ( RHO( I ) ) 4* 12*4 

12 I F ( Z ( I ) ) 13*9.13 

13 CALL SF I ELO ( D I ( I ) . T L ( I ) . A JT < I ) * R ( I ) . Z ( I ) . G ( I ) . SD * D ( I ) . FFAC . A ) 

GO TO 7 

14 I F ( NNBR ( I ) ) 201*200*201 

200 IF(I-l) 15.202.15 

202 CALL ST AT EM ( SV I *UU *D1 * D2 • P 1 . P2 • PP » PM.SS ) 

GO TO 15 

15 CALL FI NMAT (U*UU,V,MAT,A,Z( I ) ) 

DO 16 J = 1 .6 

DO 16 K = 1 * 3 

16 UU< J.K)*V( J.K) 

IF(PM) 207.1114*85 

1114 IF(FINK) 207.85*207 
85 IF(BC) 207.207.1112 
1112 IF(BRC) 207,207.208 

208 CALL BRCOR ( R3 *UU. XX.JK.VV.PM) 

MMM=N-1 

DO 209 J= 1 * 3 
DO 209 K= 1 , 3 

209 YY( J,K,MMM)=XX( J.K) 

207 GO TO 1003 

201 I F ( NNBR ( I ) -1 ) 204.204.17 

204 CALL STATEBtSVOP.W.N, YY.L.BRC .D1.D2 *PP . PM .SS .01 . 02 * 03 , F I NK ) 

17 CALL FINBRAIU.VV.V.MAT .A.Z( I) ) 

DO 18 J = 1 * 6 

DO 18 K = 1 ♦ 3 

18 VV( J.K)=V< J.K) 

I F ( NNBR ( I ) )82, 19.19 

19 IF(NNBR( I )-3) 1003,82,82 

82 CALL BRANCO ( R3 • VV * PH I ( I I ) * U ) 

N = N+ 1 
11=11+1 
MAT =0 



002800 

002810 

002820 

002830 

002840 

002850 

002860 

002870 

002880 

002890 

002900 

002910 

002920 

002930 

002940 

002950 

002960 

002970 

002980 

002990 

003000 

003010 

003020 

003030 

003040 

003050 

003060 

003070 

003080 

003090 

003100 

003110 

003120 

003130 

003140 

003150 

003160 

003170 

003180 

003190 

003200 

003210 

003220 

003230 

003240 

003250 

003260 

003270 

003280 

003290 

003300 

003310 

003320 

003330 

003340 

003350 



n n n 





BC = 1 • 

GO TO 15 


003360 




003370 


1003 


CONTINUF 


003380 

003390 




ITERATION 


003400 

003410 


2000 


CALL DELMA(UU*SVE.V.X2 .Y2.KKK. Xl *FF . PP . HH , RR . SS . Yl .REM. FINK .PM.Dl * 


003420 


1D2) 


003430 




IF(FINK) 87*87.86 


003440 


87 


IF(PM) 99*99*98 


003450 


98 


IF(KKK) 7030*7030*301 


003460 


301 


GO TO 6001 


003470 


7030 


I F ( FF- 1 *0 ) 7051 *7060*7051 


003480 


7060 


KK= 1 


003490 




M=M- 1 


003500 




GO TO 3000 


003510 


7051 


I F ( HH- 1 *0 ) 300.6001.300 


003520 


300 


IF(SS) 303.302.302 


003530 


303 


DEL ( L ) -X2-Y2, 


003540 




D1 =X 1 


003550 




D2= (X2+Y2 ) /2.0 


003560 




GO TO 61 


003570 


302 


DEL ( L ) - X 1-Y 1 


003580 




D1 = ( Xl+Yl ) /2*0 


003590 




D2 = X2 


003600 




GO TO 61 


003610 


990DEL ( L ) =V( 1 . 1 ) *V ( 2 . 2 ) *V ( 3 . 3 ) -V( 1 . 3 ) *V ( 2 . 2 ) *V ( 3 * 1 ) +V ( 1 .2 ) *V ( 2 . 3 ) *V ( 3 


003620 




1 .1 )-V( 1 *2)*V(2* 1 )*V< 3*3)+V( 1.3)*V( 3 . 2 ) * V C 2.1 )-V< 1 .1 > * V ( 3.2)*V(2 .3) 


003630 




GO TO 61 


003640 


21 


PDEL = DEL < L ) 


003650 




GO TO 25 


003660 


22 


IF(AA)23. 23*24 


003670 


23 


PDEL*DEL { L- 1 ) 


003680 




GO TO 25 


003690 


24 


PDEL r DELC 


003700 


25 


I F ( PDEL^DEL ( L ) ) 26 *26 *27 


003710 


26 


IF (AA) 28 *28 *57 


003720 


27 


BB = 0 


003730 




IF( AA) 35.35.52 


003740 


28 


AA*1.0 


003750 




BB=1.0 


003760 




DELC*DEL < L- 1 ) 


003770 




OMEGAC*OMEGA ( L-l ) 


003780 




GO TO 59 


003790 


2900MEGAM ( M ) -OMEGAC+ ( OMEGA ( L > -OMEGAC ) *ABSF < OELC ) / ( ABSF ( DELC )+ 


003800 




1 ABSF ( DEL ( L ) ) ) 


003810 




GO TO 32 


003820 


30 


0MEGA(L+1 )*OMEGA(L)-0*9*D0M 


003830 




D0M*=0* 5*D0M 


003840 


31 


L = L+ 1 


003850 




BL = L 


003860 




IF(HL-BL) 7053*7053*86 


003870 


7053 


KA= 1 


003880 




GO TO 3000 


003890 


86 


IF(JK)1000. 1000. 1002 


003900 


32 


BM=M 


003910 



83 



n n n 





I F(HM-BM) 3000*3000*33 


003920 


33 


M=M + 1 


003930 




I F ( M-3 ) 340*34*34 


003940 


34 


AD0M = ( OMEGAM ( M- 1 ) -OMEGAM ( M-2 ) ) *0 . 1 


003950 


340 


AA*0«0 


003960 




BB*0 


003970 




cc=o 


003900 




GG*0 #0 


003990 




DOM* ADOM 


004000 




DELC*0 


004010 




OMEGA ( L+l )* OMEGAM ( M-l ) +2 • *CR 


004020 




OME GAO* OMEGA ( L+ 1 ) 


004030 




GO TO 31 


004040 


35 


I F ( OMEGA ( L ) -OMEGAG ) 53 • 60 * 60 


004050 


50 


OMEGA (L+l) = OMEGA ( L ) - ( OMEGA ( L ) -OMEGAC ) / 2 . 0 


004060 




IF (OMEGA(L) -OMEGA (L+l ) -CR ) 51 *51 *31 


004070 


51 


CC=1.0 


004080 




GO TO 31 


004090 


52 


DELC=DEL ( L ) 


004100 




OMEGAC* OMEGA ( L ) 


004110 


53 


OMEGA ( L+l ) *OMEGA (L ) +DOM 


004120 




GO TO 31 


004130 


54 


IF(DELC*DEL(L) )29*29*55 


004 140 


5 500MEGAM( M ) =0MEGA ( L ) + ( OMEGA t L-l ) -OMEGA (LI) *ABSF ( DEL ( L ) I / ( ABSF < DEL ( L ) 


004150 


1)+ABSF(DEL(L-1) )) 


004160 




GO TO 32 


004170 


56 


IF ( OMEGA( L ) -OMEGAO )21*21#22 


004180 


57 


I F ( BB ) 58 # 58 * 50 


004190 


58 


BB*1.0 


004200 


59 


I F ( OMEGA ( L) -OMEGA ( L-l ) -CR ) 29*29*30 


004210 


60 


M = M-1 


004220 




GO TO 3000 


004230 


61 


IF(GG)70#70#72 


004240 


70 


I F ( DEL ( L ) ) 77 ♦ 7 1 ♦ 77 


004250 


71 


GG= 1 • 0 


004260 




OMEGA ( L+l )* OMEGA ( L ) -CR 


004270 




GO TO 31 


004280 


72 


IF(GG-1. 0)73*73*74 


004290 


73 


GG * 2 • 0 


004300 




HDEL*DEL(L) 


004310 




OMEGA ( L+l )*OMEGA(L )+2**CR 


004320 




GO TO 31 


004330 


74 


TDEL*DEL(L) *HDEL 


004340 




IF (TDEL)75*76*76 


004350 


75 


OMEGAM (M)*OMEGA(L-2» 


004360 




GO TO 33 


004370 


76 


KK= 1 


004380 




M*M- 1 


004390 




GO TO 3000 


004400 


77 


I F ( CC ) 56 * 56 * 54 


004410 






004420 




OUTPUT 


004430 

004440 


3000 


I F ( NM ) 400 #400 *403 


004450 


400 


PRINT 401 


004460 


401 


FORMAT ( 1H1// ) 


004470 



19 



PRINT 402 

402 FORMAT I15H PROGRAM VIPIPE.35X* 15HY. S. MM NHA-3.38X, 1 3H MARCH 

1 1966//.10X.108H MODAL FREQUENCIES OF A PLANAR PIPING 

2 SYSTEM ARE DETERMINED BY AN ITERATIVE PROCEDURE USING THE /29H ME 

3TH0D OF TRANSFER M ATR I CES. / / / 1 19H ******************************** 
^*************** *#***###*#****#**««***»** **#•»*****************»****. 

5******************* //) 

403 IF(NM)4058.4058.A04 

404 PRINT 405 

405 FORMAT (1H1//I 

4058 NM= 1 

4059 IFIJM4060. 4060. 4065 

4060 PRINT 4070» JJ 

40700FORMAT ( 50X ♦ 8HPR0BLEM I 2 * //35X #47HI N PLANE MODE FREQUENCIES (RADIA 
INS PER SECOND) / /45X .4HM0DE # 1 5X # 10H FREQUENCY /) 

GO TO 4079 
4065 PRINT 4075 » JJ 

40750FORMATI 50 X * 8 HP ROB LEM I 2 # / / 33X # 5 1H0UT OF PLANE MODE FREQUENCIES (RA 
1DIANS PER SECOND) //45X *4HM0DE » 15X • 10H FREQUENCY /) 

4079 MH=M 

PRINT 408 # ( M*0MEGAM( M ) #M* 1 #MH ) 

408 FORMAT ( 46X » I 2 ♦ 1 7X # F 1 4 . 8 /) 

IF(KA) 5000 #5000# 5010 

5010 PRINT 5020 

5020 FORMAT ( 10X * 53HPROBLEM TERMINATED DUE TO OMEGA STORAGE LIMITATION. 

1 /) 

5000 IF(KK)4080*4090*4080 

4080 PRINT 4085 

40850FORMAT( 10X * 69HPR0BLEM TERMINATED DUE TO SIGNIFICANT FIGURE LIMITAT 
1 I ON OF COMPUTER. / )’ 

4 09 0 IF (AMYK) 41 2 #409*412 

409 PRINT 411 

4 1 10F0RMAT (16X*53H A LUMPED MASS APPROACH WAS EMPLOYED# A 

1ND ) 

GO TO 414 

412 PRINT 413 

413 FORMAT (24X*45HA DISTRIBUTED MASS APPROACH WAS EMPLOYED AND ) 

414 IF(SD)415*415#420 

415 I F ( R I ) 4 16*416*418 

416 PRINT 417 

4 1 70FORMAT (2X#67HTHE EFFECTS OF SHEAR DEFLECTION AND ROTARY INERTIA W 
1 ERE CONSIDERED. /) 

GO TO 4610 
418 PRINT 419 

41 90FORMAT (3X#63HTHE EFFECTS OF SHEAR DEFLECTION WERE CONSIDERED WHIL 
IE THOSE OF /3X#35HR0TATI0NAL INERTIA WERE NEGLECTED. /) 

GO TO 4610 

420 IF(RI )421 #421 *423 

421 PRINT 422 

4220F0RMAT (3X»72HTHE EFFECTS OF ROTATIONAL INERTIA WERE CONSIDERED WH 
1ILE THOSE OF SHEAR /3X*27HDEFLECTION WERE NEGLECTED. /) 

GO TO 4610 
423 PRINT 424 

4240F0RMATI 3X*66HTHE EFFECTS OF SHEAR DEFLECTION AND ROTARY INERTIA WE 
IRE NEGLECTED. /) 

4610 IF(PM) 425*425»4713 



004480 

004490 

004500 

004510 

004520 

004530 

004540 

004550 

004560 

004570 

004580 

004590 

004600 

004610 

004620 

004630 

004640 

004650 

004660 

004670 

004680 

D04690 

004700 

004710 

004720 

004730 

004740 

004750 

004760 

004770 

004780 

004790 

004800 

004810 

004820 

004830 

004840 

004850 

004860 

004870 

004880 

004890 

004900 

004910 

004920 

004930 

004940 

004950 

004960 

00A970 

004980 

004990 

005000 

005010 

005020 

005030 



90 



4713 IF(SS) 4611.4612,4612 
4611 I F ( PP- 1 • ) 4613,4614,4615 



4613 


PRINT 4623 










4623 


FORMAT ( 24X , 35HPESTEL 
GO TO 425 


MAHRENHOLTZ 


METHOD! A) 


USED. 


/) 


4614 


PRINT 4624 










4624 


FORMAT! 24X,35HPESTEL 
GO TO 425 


mahrenholtz 


METHOD! B) 


USED. 


/) 


4615 


PRINT 4625 










4625 


FORMAT! 24X,35HPESTEL 
GO TO 425 


MAHRENHOLTZ 


METHOD ( C ) 


USED. 


/) 


4612 


I F ( SS- 1 • ) 4616,4617, 


4618 








4616 


PRINT 4626 










4626 


FORMAT! 24X,35HPESTEL 
GO TO 425 


MAHRENHOLTZ 


METHOD ( D ) 


USED.- 


/) 


4617 


PRINT 4627 










4627 


FORMAT! 24X,35HPESTEL 
GO TO 425 


MAHRENHOLTZ 


METHOD! E) 


USED. 


/) 


4618 


PRINT 4628 










4628 


FORMAT! 24X,35HPESTEL 


MAHRENHOLTZ 


METHOD! F) 


USED. 


/) 


425 


IF ( PM) 9010*9020t9000 








9010 


PRINT 9012 










9012 


FORMAT! 24X.20HDELTA 
GO TO 9000 


METHOD USED. 


/ ) 






9020 


PRINT 9022 










9022 


FORMAT! 24X.29HM0DI FI 


ED DELTA METHOD USED. 


/) 




9000 


PRINT 440 











4400 FORMAT ( 1 18H ****************************************************** 

1 ************************************************** ******* ***** j 

I F ( GDABC- 1 • ) 45 11, 452 Ot 4401 
4401 PRINT 4402 

44020FORMAT ( 35H SECTION LENGTH AND SUBSECTION DATA/23X , 14HSECT ION NUMBE 
1R.8X.25HLENGTH OF SECT I ON ( INCHES ) , 5X . 2 1HNUMBER OF SECTIONS ) 

PRINT 443, ( I »TL( I ) ,Z( I ) ,1*1, NS ) 

443 FORMAT ( 29X , 1*3, 23X . F7.2 ,2lX ,F4.0 ) 

LM = L 

IF! PM) 7054,7054,7055 
7055 PRINT 441 

44 10F0RMAT (11H GRAPH DATA / 22 X ♦ 16H I TERAT I ON NUMBER * 4X .30HFREQUENCY (R 
1ADIANS PER SECOND) ,3X,30H VALUE OF REMAINDER ) 

GO TO 7057 
7054 PRINT 7056 

70560FORMAT ( 1 1H GRAPH DATA /22X , 16HI TERATION NUMBER . 4X , 30HFRE0UENCY (R 
1 ADI ANS PER SECOND) *3X,30HVALUE OF FREQUENCY DETERMINANT ) 

7057 PRINT 444, ( L .OMEGA (L) , DEL ( L ) ,L*1 ,LM) 

444 FORMAT ( 29X , I 3*23X ,F1 1 .6, 18X , E15. 8 ) 

4520 PRINT 442 

442 FORMAT ( IX, 1 OH IN PUT DATA /) 

PRINT 4521 

452 10FORMAT ( IX , 9HC0MP0NENT ,4X ,8HD I AMETER , 5X , 14HWALL THICKNESS ,5X,6HLEN 
1GTH,5X,9HRADIUS OF ,9X , 8H I NCLUDED ♦ 5X , 7HDENS ITY , 5X , 7HELAST IC . 5X , 5HSH 
2EAR /57X,9HCURVATURE,5X,12HANGLE OF ARC • 16X ♦ 7HM0DULUS , 5X »7H MODULUS 
3 /) 

DO 4525 1 = 1, NS 
I F ( P ( I ) )4524, 4524, 4523 

452 4 PRINT 4 522,1 ,D< 1 ) *DT ( I ) , TL f I ) ,RHO( I ) , THETR ( I ) ,BAMU( I ) ,E( I ) ,G( I ) 



005040 

005050 

005060 

005070 

005080 

005090 

005100 

005110 

005120 

005130 

005140 

005150 

005160 

005170 

005180 

005190 

005200 

005210 

005220 

005230 

005240 

005250 

005260 

005270 

005280 

005290 

005300 

005310 

005320 

005330 

005340 

005350 

005360 

005370 

005380 

005390 

005400 

005410 

005420 

005430 

005440 

005450 

00546.0 

005470 

005480 

005490 

005500 

005510 

005520 

005530 

005540 

005550 

005560 

005570 

005580 

005590 



91 



4 5220FORMAT !3X.I3.8X.F7.3.10X.F6.3.7X.F8.3.4X.F8.3*9X.F6.2.8X*F7.2.4X. 
1E8.2.4X.E8.2 ) 

GO TO 4525 
4523 SH = 1 • 

4525 CONTINUE 

IF ( SH ) 4530 #4530 #4526 

4526 PRINT 4532 

4532 FORMAT ( / / IX . 7 HH ANGERS /) 

PRINT 4527 

4 52 70FORMAT ( 2X * 9HCOMPONENT * 14X . 3HCLX . 14X . 3HCLY . 1 4X . 3HCLZ • 14X . 3HCTX . 14X 
1 • 3HCT Y . 14X.3HCTZ /) 

4528 DO 4530 1*1. NS 

I F ( P ( I ) >4530.4530*4529 

4529 PRINT 453 1 . I . CLX ( I > . CLY Cl) . CLZ ( I > . CTX ( I) . CTY ( I ) . CTZ ( I > 

4 531 FORMAT! 3X . I 3 V 1 5X . F10 . 2 . 7X . F 10 . 2 . 7X . F 10 . 2 • 7X. F10 • 2 . 7 X • F10 • 2 . 7X • F 10 . 
12 ) 

4530 CONTINUE 
PRINT 448 

448 FORMAT L//20H BOUNDARY CONDITIONS > 

PRINT 449. (SVI ( I .1 > . 1*1 .6) 

449 FORMAT (5H SVI=6F3.0) 

PRINT 451 . ( SVE ( I .1 >. 1 = 1.6) 

451 FORMAT ( 5H SVE=6F3*0) 

IF(NBR)4511. 4511. 4502 
4502 IF! JK)4505. 4505. 4515 

4505 PRINT 4510. ( (SVIP! I . J) . 1=1 .6) . J=2.NNN> 

4510 FORMAT (6F3.0) 

GO TO 4511 

4515 PRINT 4517. ( ( SVOP ( I *J) *1=1 .6) .J=2.NNN) 

4517 FORMAT (6F3.0) 

4511 CONTINUE 

C SORT 

I F ( GRAPH- 1 • ) 6001.6000.6001 
6000 LM1*LM-1 

DO 6100 L* 1 * LM 
EDL ( L ) *DEL( L ) 

6100 CONTINUE 

DO 6600 L»1.LM1 
LP1 *L+ 1 

DO 6600 M«LP1.LM 

IF ( OMEGA! L) “OMEGA! M) ) 6600.6600.6003 
6003 TEMP=OMEGA( L) 

OMEGA!L)=OMEGA(M) 

OMEGA ( M ) =TEMP 
TEMP=EDL ( L ) 

EDL ( L ) =EDL ( M ) 

EDL ( M ) *TEMP 
6600 CONTINUE 
C GRAPH 

NUMPTS*60 
K = 0 

DO 6006 J= 1 . 12 
6006 I T I TLE ( J ) =8H 

I T I TLE ( 3 ) =8H KIM Y S 
I T I TLE ( 8 ) =8HFREQUENC 
I T I TLE ( 9 ) =8HY VS R 



005600 

005610 

005620 

005630 

005640 

005650 

005660 

005670 

005680 

005690 

005700 

005710 

005720 

005730 

005740 

005750 

005760 

005770 

005780 

005790 

005800 

005810 

005820 

005830 

005840 

005850 

005860 

005870 

005880 

005890 

005900 

005910 

005920 

005930 

005940 

005950 

005960 

005970 

005980 

005990 

006000 

006010 

006020 

006030 

006040 

006050 

006060 

006070 

006080 

006090 

006100 

006110 

006120 

006130 

006140 

006150 



92 



I r I TLE ( 10) =8HEMA1NDER 
LTl=LM/NUMPTS 
IF(LTl) 6024,6025.6024 
6025 NUMPTS=LM 
LABEL=4H 



CALL DR AW ( NUMPTS. OMEGA. EDL. 0.0. LABEL. ITITLE.0,0.0.0. 
110, K. LAST) 

GO TO 6001 



0.0.9 . 



602 A DO 6020 J=1.17 • 

LT2*0 

LT2*LT2+J*NUMPTS 
I F ( LT2-LM ) 6020.6021.6022 

6020 CONTINUE 

6021 LT=LT1 

GO TO 6023 

6022 LT=LTl+ 1 
numend=lm-lti*numpts 



6023 L=1 

DO 7021 1 = 1, LT 
I F ( I — 1 ) 600 8,6041,6008 
6008 IF(I-LT) 8011,6011.6011 

8011 I F ( I -3 ) 6042.6043,8012 

8012 IFIl-5) 6044.6045.8013 

8013 I F < I —7 ) 6046,6047,8014 

8014 I F ( I -9 ) 6048.6049.8015 

8015 IF(I-LT) 6050,6011,6001 

6041 LABEL=4H ONE 
GO TO 6009 

6042 LABEL=4H TWO 
GO TO 6010 

6043 LABEL=4HTHRE 
GO TO 6010 



6044 LABEL=4HFOUR . 
GO TO 6010 

6045 LABEL=4HF I VE 
GO TO 6010 

6046 LABEL=4H SIX 
GO TO 6010 

6047 LABEL=4HSEVN 
GO TO 6010 

6048 LABEL=4H EIT 
GO TO 6010 

6049 LABEL*4HNINE 



GO TO 6010 
6050 LABEL=4H TEN 
GO TO 6010 



6009 

6031 

6027 



0 6031 N* 1 .NUMPTS 
DL ( N) =EDL( L+N-l ) 

MEGA ( N ) =OMEGA ( L+N- 1 ) 

ORMAT^lHl »///// / / SXWHNUMPTSs ,I3,5X,7HEDL(1) = »E12«9) 

:ALL DRAW (NUMPTS, OMEGA. EDL, 0.0. LABEL, ITITLE.O. 0.0. 0.0.0. 



9. 



110, K. LAST ) 
GO TO 7021 
6010 L“L+NUMPTS 



006160 

006170 

006180 

006190 

006200 

006.210 

006220 

006230 

006240 

006250 

006260 

006270 

006280 

006290 

006300 

006310 

006320 

006330 

006340 

006350 

006360 

006370 

006380 

006390 

006400 

006410 

006420 

006430 

006440 

006450 

006460 

006470 

006480 

006490 

006500 

006510 

006520 

006530 

006540 

006550 

006560 

006570 

006580 

006590 

006600 

006610 

006620 

006630 

006640 

006650 

006660 

006670 

006680 

006690 

006700 

006710 





OMEGAX=OMEGA ( L ) 


006720 




DO 6030 N=1.NUMPTS 


006730 




EDL ( N ) =EDL ( N+L-l ) 


006740 




OMEGA ( N ) “OMEGA ( L+N- 1 ) 


006750 


6030 


OMEGA) N ) “OMEGA! N) -OMEGAX 


006760 




PRINT 6032.NUMPTS.EDL( 1 ) 


006770 


6032 


FORMAT < ///////5X.7HNUMPTS* . 13 , 5X , 7HEDL < 1 ) = , E 1 2 .9 ) 


006780 




CALL DRAW (NUMPTS. OMEGA . EDL .0 . 0 . LABEL . I T I TLE . 0 .0 .0 .0 .0*0.9. 


006790 




110. K. LAST) 


006800 




GO TO 7021 


006810 


6011 


L=L+NUMPTS 


006820 




OMEGAX“OMEGA ( L ) 


006830 




DO 6033 N»1,NUMEND 


006840 




EDL ( N) =EDL ( N+L-l ) 


006850 




OMEGA (N (“OMEGA t N+L- 1 ) 


006860 


6033 


OMEGA ( N ) “OMESA ( N ) -OMEGAX 


006870 




PRINT 6051.NUMEND.EDLI1) 


006880 


6051 


FORMATI///////5X. 7HNUMEND= .13. 5X.7HEDLI 1 > “.E12.91 


006890 




LABEL* AHLAST 


006900 




CALL DR AW (NUMEND. OMEGA. EDL, 0.0. LABEL, I T I TLE . 0 ,0 .0 .0 . 0 .0 . 9 . 


006910 


1 10 t K * LAST ) 


006920 




GO TO 6001 


006930 


7021 


CONTINUE 


006940 


6001 


IF(HH-l.O) 8024*8023 *8023 


006950 


8023 


MH = M 


006960 




LM e L 


006970 


8024 


DO 446 M* 1 * MH 


006980 


446 


OMEGAMI M) *0 *0 


006990 




DO 447 L a 1 • LM 


007000 




OMEGA < L ) = 0# 0 


007010 




EDL ( L ) *0* 


007020 


447 


DEL ( L ) *0* 0 


007030 




DELC=0.0 


007040 




OMEGAC a O.O 


007050 




IF(HH-l.O) 8001*8022*8001 


007060 


8001 


IF(PM) 8031.8031*8033 


007070 


8033 


IF(REM) 8031.8032.8032 


007080 


8031 


IF (MN) 8022*8022 *452 


007090 


8032 


I F ( REM- 1 • ) 8041*8042*8043 


007100 


8041 


IF(PP-l) 8043*8043*8072 


007110 


8072 


IF(SS) 8073*8074.8074 


007120 


8073 


ss*o. 


007130 




GO TO 8022 


007140 


8074 


IF(SS-1* ) 8042. 8042. 8075 


007150 


8075 


SS r -l* 


007160 




PP*0. 


007170 




REM*-1# 


007180 




GO TO 8033 


007190 


8042 


SS=SS-1* 


007200 




IF(SS) 8051 * 8052 ♦ 805 3 


007210 


80 51 


SS«1. 


007220 




GO TO 8022 


007230 


8052 


SS = 2 • 


007240 




GO TO 8022 • 


007250 


8053 


SS=-1* 


007260 




REM=-1. 


007270 





GO TO 8033 


007280 


8043 


PP®PP- 1. 


007290 




IF(PP) 8061 .8062.8063 


007300 


8061 


PP = 1* 


007310 




GO TO 8022 


007320 


8062 


PP = 2 • 


007330 




GO TO 8022 


007340 


8063 


PP®0* 


007350 




REM®-1 * 


007360 




GO TO 8033 


007370 


452 


PP® 0. 


007380 




RR* 0 * 


007390 




SS=-1. 


007400 


8022 


L® 1 


007410 




M® 1 


007420 




AA*0.0 


007430 




BB®0.0 


007440 




CC=0*0 


007450 




FF® 0 *0 


007460 




GG=0.0 


007470 




KK®0 


007480 




KA®0 


007490 




KKK = 0 


007500 




ADOM®BDOM 


007510 




CR = BR 


007520 




OMEGA( 1 )®OMEGAI 


007530 




OMEGAO«OMEGAI 


007540 




DOM®ADOM 


007550 




IF(HH-l.O) 8021 *7040 * 802 1 


007560 


7040 


IF(SS) 8191*8192 #8 192 


007570 


8191 


I F ( PP- 1 • ) 8091.8092*8093 


007580 


8192 


I F ( SS- 1 • ) 8291.8292*8293 


007590 


8091 


PRINT 8081 


007600 


8081 


FORMAT (1H1.////18HMETHOD A FAILED* ///) 


007610 




GO TO 8021 


007620 


8092 


PRINT 8082 


007630 


8082 


FORMAT ( 1H1.////18HMETHOD B FAILED* ///) 


007640 




GO TO 8021 


007650 


8093 


PRINT 8084 


007660 


8084 


FORMAT ( 1H1 • ////l 8HMETHOD C FAILED* ///) 


007670 




GO TO 8021 


007680 


8291 


PRINT 8281 


007690 


8281 


F0RMAT(1H1.////18HMETH0D D FAILED* ///) 


007700 




GO TO 8021 


007710 


8292 


PRINT 8282 


007720 


8282 


FORMAT ( 1H1.////18H METHOD E FAILED* ///) 


007730 




GO TO 8021 


007740 


8293 


PRINT 8283 


007750 


8283 


FORMAT ( 1H1*////18HMETH0D F FAILED* ///) 


007760 


8021 


IF(PM) 453.453.8035 


007770 


8035 


IF(REM) 8037*8036.8036 


007780 


8036 


IF(JK) 998*999.998 


007790 


8037 


REM=ERM 


007800 


453 


IF(PM) 6081.6082.6080 


007810 


6081 


IF(PMM-1.) 6084.6080.6084 


007820 


6084 


I F ( PMM*0 • 1-10*1) 6086.6085.6086 


007830 



KJ KJ 



6085 PM* 1 • 

GO TO 8036 

6086 PM=0. 

6082 ?F (PMM*0. 1-1. 1 ) 6080.6080.6085 
6080 I F t MN ) 9001 .9001 .8034 
8034 MN*0 

PMxPPM 
GO TO 998 
9001 JJ=JJ+1 

DO 454 1*1. NS 
P( I) *0. 

454 Z< 1 )=0.0 
GO TO 450 
600 END 

SUBROUTINE SUBSEC I RHO . THE TA , D , TL . HM.Z .OMEGA I ,HMH ) 
RHOD*ABSF(RHO) 

IF(HMH) 41*40*41 



007840 

007850 

007860 

007870 

007880 

007890 

007900 

007910 

007920 

007930 

007940 

007950 

007960 

007970 

007980 

007990 

008000 

008010 

008020 

008030 

008040 



40 I F ( HM-5 .0 18.8.7 


008050 


8 I F ( OMEGA I -800. ) 10*T.7 


008060 


7 HM=6. 


008070 


GO TO 10 


008080 


41 HM=HMH 


008090 


10 RA=TL/D 


008100 


IFIRHOD-0. >3.3.4 


008110 


3 IF ( RA-1. ) 5. 5.6 


008120 


5 Z = 0. 


008130 


RFTURN 


008140 


6 IFIRA-3.) 15.15.16 


008150 


15 Z = 2.0 


008160 


RETURN 


008170 


16 IFIRA-6. ) 17. 17. 18 


008180 


17 Z*3.0 


008190 


RETURN 


008200 


18 IFCHM-3.il. 1.2 


008210 


1 Z=6.0 


008220 


RETURN 


008230 


2 Z*2.0»HM 


008240 


RETURN 


008250 


4 rd=rhod/d 


008260 


IFCRO-1.) 19.19.23 


008270 


19 Z*0 


008280 


RETURN 


008290 


23 IF(RD-3. >24,24,25 


008300 


24 IFCTHET A- .44)26.26. 27 


008310 


26 Z*0. 


008320 


RETURN 


008330 


27 IFCTHET A- 1.57 >28.28. 29 


008340 


28 Z*4.0 


008350 


RETURN 


008360 


29 Z=6.0 


008370 


RETURN 


008380 


25 IFCRD-6. >30,30.31 
30 IF ( THETA-. 44) 32 .32 .33 


008390 



4 



96 



n n 



32 


2 = 4.0 


008400 




RETURN 


008410 


33 


IF( THETA-1. 57)34.34, 35 


008420 


34 


2=6.0 


008430 




RF TURN 


008440 


35 


2 = 8.0 


008450 




RETURN 


008460 


31 


IF ( THETA- .44 ) 36.36 ,3 


008470 


36 


2 = 4.0 


008480 




RETURN 


008490 




END 


008500 






008510 

008520 




SUBROUTINE D I STM ( AMU , A I Y , OMEGA ,TL,R,D,DI , G, SD , RI ,E,FFAC,U) 


008530 




DIMENSION U ( 6 « 6 ) 


008540 




ARA= ( 3* 141592654/4.0)* ( D**2-DI **2 ) 


008550 




AREA=ARA*FFAC 


008560 




BE=TL*OMEGA*SQRTF< AMU/ (E*ARA) ) 


008570 




P4=R*AMU*TL**4*OMEGA**2 


008580 




IF(RI)1,1,2 


008590 


1 


T=R*AMU* ( A I Y*OMEGA*TL ) **2 


008600 




GO TO 3 


008610 


2 


T = 0 • 0 


008620 


3 


I F ( SD ) 4 ,4 , 5 


008630 


4 


S= ( AMU* ( OMEGA *TL ) **2 ) / ( G*AREA ) 


008640 




GO TO 6 


008650 


5 


s=o.o 


008660 


6 


SPT = S+T 


008670 




VV=SORTF( P4+(S-T) **2/4.0) 


008680 




AL1=SQRTF(VV-.5*SPT) 


008690 




AL2=SGRTF(VV+.5*SPT) 


008700 




CL=1.0/(AL1**2+AL2**2) 


008710 




CHL 1 = ( EXPF ( 4L1)+1./EXPF(AL1) )/2. 


008720 




SHL 1 = ( EXPF ( AL 1 ) -1 • /EXPF t AL 1 ) )/?. 


008730 




CL2 = COSF( AL2 ) 


008740 




SL2 = SI NF ( AL2 ) 


008750 




SBE = SI NF ( BE ) * 


008760 




CBE =COSF ( BE ) 


008770 




COO=CL* (AL2**2*CHL1+AL1**2*CL2 ) 


008780 




C01=CL* ( ( AL2**2*SHL1) / AL 1 + ( AL 1 **2*SL2 ) / AL2 ) 


008790 




C02 = CL* (CHL1-CL2) 


008800 




C03 =CL* (SHL1/AL1-SL2/AL2) 


008810 




DO 7 J= 1 , 6 


008820 




DO 7 K= 1 , 6 


008830 


7 


U( J,K)=0.0 


008840 




U ( 1 , 1 ) =CBE 


008850 




Ut 1,6)=TL*SBE/(BE*ARA*E) 


008860 




U(2,2)=CO0~S*CO2 


008870 




U( 2,3)=TL*(C01-SPT*C03) 


008880 




U(2,4)=R*C02*TL**2 


008890 




U(2,5)=(R*TL**3/P4)*<( P4+S**2 ) *C03-S*C0 1 ) 


008900 




U( 3,2) =P4*C03/TL 


008910 




U( 3,3) =CO0-T*CO2 


008920 




U ( 3 ,4 ) =TL*R* ( C01-T*C03 ) 


008930 




U( 3,5)=U( 2,4) 


008940 




RR=1.0/R 


008950 



9T 



n n 



2 

3 

4 

5 

6 
7 



1 



U< 4*2) =P4*RR*C02/TL**2 
U(4»3)=RR*( ( P4 + T ** 2 ) *C03-T*C01 ) /TL 
U ( 4 * 4 ) =U ( 3*3) 

U ( 4 * 5 ) =U ( 2*3) 

U( 5*2)=P4*RR* (C01-S*C03 )/TL**3 
U( 5 ♦ 3 ) =U (4*2) 

U( 5 *4 ) =U ( 3 * 2 ) 

U(5*5)=U(2*2) 

U( 6* T) *-AMU*TL*0MEG4**2*SBE/BE 

U ( 6 * 6 ) = CBE 

RETURN 

END 



SUBROUTINE DISTMO ( AMU * A I X . A I Y * OMEGA * T L * A JT * R * D » D 1 *G»SD»RI *FFAC*U) 
DIMENSION U ( 6 *6 ) 

W= 1 • / ( G*A JT ) 

AREA= ( 3.141592654/4,0 ) *<D**2-DI **2 )*FFAC 
BE=+SQRTF ( AMU * < T L *A I X*OMEGA ) **2*W ) 

P4=R*AMU * TL **4*OMEGA**2 
IF( RI ) 2.2*3 

T =R*AMU* ( A I Y*OMEGA*TL ) **2 
GO TO 4 
T = 0 • 0 

IF(SD) 5*5*6 

S=( AMU* ( OMEGA *TL)**2 ) / ( G*AREA ) 

GO TO 7 
S*0.0 

VV = + SGRTF (P4+< S-T )** 2/4.0 ) 

SPT =S+T 

ALl*+SGRTF< VV-.5*SPT ) 

AL2 c +SORTF( VV+. 5* SPT ) 

CL*1«0/(AL1**2+AL2**2) 

CHL1*(EXPF(AL1 )+l*/EXPF (AL1 ) ) / 2 • 

SHL1*( EXPF ( AL1 ) -1 • /EXPF ( AL 1 ) )/2. 

CL2 *COSF ( AL2 ) 

SL2*SINF( AL2 ) * 

SBE*SI NF ( BE ) • 

CBE*COSF ( BE ) 

COO=CL*(AL2**2*CHL1 + AL1**2*CL2) 

C01*CL*< ( AL2**2*SHL1 )/ALl + ( AL 1 **2*SL2 ) / AL2 ) 

C02=CL*(CHL1-CL2 ) 

C03=CL*<SHL1/AL1-SL2/AL2 ) 

DO 1 J= 1 • 6 

DO 1 K*l#6 

U( JtK) *0.0 

U ( 1 * 1 ) =CBE 

Ul 1 t2)«TL*W*SBE/BE 

U(2*1)=-AMU * TL *AIX**2*S8E/BE 

U ( 2 . 2 ) *CBE 

U( 3*3) s C00“S # C02 

U ( 3 * 4 ) * TL * ( CO 1-SPT *C03 ) 

U( 3*5) =R*C02*TL**2 

U( 3.6) »<R*TL**3/P4)*(-S*C01+< P4+$**2 )*C03) 

U(4*3) *P4*C03/TL 
U(4*4)=COO-T*C02 



008960 

008970 

008980 

008990 

009000 

009010 

009020 

009030 

009040 

009050 

009060 

009070 

009080 

009090 

009100 

009110 

009120 

009130 

009140 

009150 

009160 

009170 

009180 

009190 

009200 

009210 

009220 

009230 

009240 

009250 

009260 

009270 

009280 

009290 

009300 

009310 

009320 

009330 

009340 

009350 

009360 

009370 

009380 

009390 

009400 

009410 

009420 

009430 

009440 

009450 

009460 

009470 

009480 

009490 

009500 

009510 



U ( 4 * 5 ) = TL*R* (C01-T*C03 ) 


009520 


Ut 4,6) =U( 3.5) 


009530 


RR=1.0/R 


009540 


U(5.3)=P4*RR*C02/TL**2 


009550 


U(5.4)=RR*(-T*C01+(P4+T**2)*C03)/TL 


009560 


U(5»5)=U(4.4) 


009570 


Ut5.6)=Ut3,4) 


009580 


Ut 6.3)=P4*RR* (C01-S*C03)/TL**3 


009590 


U(6»4)=U(5»3) 


009600 


U(6»5)=U(4»3) 


009610 


U(6,6)=U(3,3) 


009620 


RETURN 


009630 


END 


009640 




009650 




009660 


SUBROUTINE RIGID ( AMU »TL , OMEGA .AIY.RI.U) 


009670 


DIMENSION U ( 6 * 6 ) 


009680 


AM=AMU*TL 


009690 


DO 1 J=1 »6 


009700 


DO 1 K=1 *6 


009710 


U( J»K)=0.0 


009720 


U( 1.1)=1.0 


009730 


U( 6 * 1 ) =-AM*OMEGA**2 


009740 


U<6. 61=1.0 


009750 


U<2»2>=1.0 ^ 


009760 


U( 2 * 3 ) =TL 


009770 


U( 3.31*1.0 


009780 


U(4»2) =AM*TL»OMEGA**2/2.0 


009790 


IFfRI >2*2.3 


009800 


U(4.3)=AM*0MEGA**2*( TL*»2/6.0-A IY**2 ) 


009810 


GO TO 4 


009820 


U( 4.3) =AM*(OMEGA*TL ) **2/6.0 


009830 


U t 4*4) »1 « 0 


009840 


U( 4 * 5 ) =TL 


009850 


U( 5 * 2 ) =-U (6*1) 


009860 


U( 5*3) =U( 4,2) 


009870 


U(5.5)=1.0 


009880 


RETURN 


009890 


END 


009900 

009910 

009920 


SUBROUTINE RIGIO ( AMU * TL. A I X .OMEGA * A I Y * AJT *SDR I *U ) 


009930 


DIMENSION UI6.6) 


00994G 


DO 2 J= 1,6 


009950 


DO 2 K=1 *6 


009960 


Ut J*K)«0.0 


009970 


AM=AMU*TL 


009980 


Utl,l)=1.0 


009990 


U( 2 . 1 ) »-AMU*TL* ( A I X*OMEGA ) **2 


010000 


Ul 2.2) *1.0 


010010 


U(3,3)=1.0 


010020 


U( 3 » 4) =TL 


010030 


U(4,4)*1.0 


010040 


U ( 5 . 3 ) =AM*TL*OMEGA**2/2 .0 


010050 


IFtSDRI-O.) 30.30.31 


010060 


Ut 5 .4) =AM*0MEGA**2*t TL**2/6.0-AIY**2 ) 


010070 



99 



on 



GO TO 1 

31 U( 5.4) = AM*OMEGA** 2 * ( TL**2/6*0) 
1 U( 5 • 5 ) = ] .0 
U ( 5 . 6 ) = TL 

U ( 6 * 3 ) =AM*0MFGA**2 
U(6.4)=U<5.3) 

U( 6 1 6 ) = 1 *0 

RETURN 

END 



SUBROUTINE STiFCO ( THE TA * RHO .U ) 
DIMENSION U ( 6 * 6 ) 

DO 1 J= 1 *6 
DO 1 K'1.6 

1 UUtK)=0.0 

I F ( RHO ) 2 » 2 . 3 

2 V=-1*0 
GO TO 4 

3 V=1.0 

4 CT s COSF (THETA ) 

ST = $INF ( THETA ) 

U( 1*1) =CT 

U( 1.2)»-V*ST 
U ( 2 * 1 ) * V*ST 
U ( 2 • 2 ) =CT 
U( 3*3 ) *1*0 
U(4#4)=l«0 
U(5.5)»CT 
U ( 5 . 6 ) s -V*ST 
U ( 6 * 5 ) = V*ST 
U ( 6 *6 ) =CT 
RETURN 
END 



SUBROUTINE STIFOO ( THETA .RHO .U ) 
DIMENSION U ( 6 * 6 ) 

IF (RHO) 1*1*2 

1 V=-1.0 
GO TO 3 

2 V=1.0 

3 CT = COSF ( THETA ) 

ST C S I NF ( THE T A ) 

DO 4 J- 1 . 6 

DO 4 K = 1 . 6 

4 U( J .K) =0*0 
U C 1 * 1 ) =CT 
U( 1.4)*ST*V 
U ( 2 • 2 ) =CT 

U ( 2 . 5 ) =ST*V 
U(3.3)=l*0 
U ( 4 . 1 ) »-ST*V 
U ( 4 . 4 ) = CT 
U( 5.2) =-ST*V 
U( 5.5) =CT 



010080 
010090 
010100 
010110 
010120 
010130 
010 140 
010 150 
010160 
010170 
010180 
010 190 
010200 
010210 
010220 
010230 
010240 
010250 
010260 
010270 
010280 
010290 
010300 
010310 
010320 
010330 
010340 
010350 
010360 
010370 
010380 
010390 
010400 
010410 
010420 
010430 
010440 
010450 
010460 
010470 
010480 
010490 
010500 
010510 
010520 
010530 
010540 
010550 
010560 
010570 
010580 
010590 
010600 
010610 
010620 
010630 






i 



n n n n 



U(6*6)=1.0 

RETURN 

END 



SUBROUTINE ST AVEC < SV . N * BC ) 
DIMENSION BC ( 6 * 22 ) 

DO 10 J=1 • 6 
10 BC(J»N)=0. 

I F ( SV- 1 • ) 1 1 1 * 2 

1 BC ( 4» N ) = 1 • 0 
BC(5tN)-l*0 
BC(6tN)=1.0 
RETURN 

2 I F ( SV-2 • ) 3 * 3 * 4 

3 BC(ltN)=1.0 
BO 2 tN) =1.0 
BC ( 3 * N ) * 1 • 0 
RETURN 

4 I F ( SV- 3 • ) 5 * 5 # 6 

5 BC(3»N)=1.0 
BC(5tN)=1.0 
BC(6*N)=1.0 
RETURN 

6 I F ( SV— 4 • ) 7 • 7 * 8 

7 BC(ltN)=1.0 
BC ( 3 • N ) = 1 • 0 
BC ( 5 »N ) = 1 • 0 
RETURN 

8 BC(2*N)=1.0 
BC(3tN)=1.0 
BC(6#N)=1.0 
RETURN 

END 



SUBROUTINE STAVED (SV*N*BC) 
DIMENSION BC ( 6 t 22 ) 

DO 10 L=1 »6 
10 BC ( L # N ) =0 . 0 

I F ( S V— 1 • 0 ) 1 * 1 * 2 

1 BC(2#N)=1.0 
BC ( 5 * N ) *1 • 0 
BC(6*N)*1.0 
RETURN 

2 I F ( SV-2 . 0 ) 3 * 3 *4 

3 BCU»N)=1.0 
BC(3»N)M.0 
BC(4tN) =1.0 
RETURN 

4 I F ( SV-3 .0 ) 5 t 5 *6 

5 BC(2»N)»1.0 
BC(4tN)=1.0 
BC(6tN)*1.0 
RETURN 

6 I F ( SV-4 .0 ) 7 * 7 ♦ 5 



010640 

010650 

010660 

010670 

010680 

010690 

010700 

010710 

010720 

010730 

010740 

010750 

010760 

010770 

010780 

010790 

010800 

010810 

010820 

010830 

010840 

010850 

010860 

010870 

010880 

010890 

010900 

010910 

010920 

010930 

010940 

010950 

010960 

010970 

010980 

010990 

011000 

011010 

011020 

011030 

011040 

011050 

011060 

011070 

011080 

011090 

011100 

011110 

011120 

011130 

011140 

011150 

011160 

011170 

011180 

011190 



101 



no r> n 



7 


BC( 1 .Nl =1.0 


011200 




BC(4.N)=1.0 


011210 




BC(6.N)=1.0 


011220 




RETURN 


011230 




END 


Oil 240 






011250 

011260 




SUBROUTINE INVERT (A.N.D.L.M) 


011270 




PROGRAM FOR FINDING THE INVERSE OF A NXN MATRIX 


011280 






011290 




DIMENSION A ( 25.25 ) .L ( 25 ) ,M(25) 


011300 




SEARCH FOR LARGEST ELEMENT 


011310 




D=1 .0 


011320 




D080 K = 1 »N 


011330 




L ( K ) =K 


011340 




M(K)=K 


011350 




8 1 GA = A ( IC» K ) 


011360 




D020 I=X.N 


011370 




D020 J = X »N 


011380 




IF(ABSF(BIGA)-ABSF(A( I ,J) ) ) 10.20.20 


011390 


10 


BIGA=A ( I .J) 


011400 




L I X ). = I 


011410 




M( X)*J 


011420 


20 


CONTINUE 


011430 




INTERCHANGE ROWS 


0 1 1440 




J=L ( X 1 


011450 




IF ( L ( X ) -X ) 35.35.25 


011460 


25 


D030 1=1. N 


011470 




HOLD=-A ( X » I ) 


011480 




A ( K. • I )»A( J.I ) . 


011490 


30 


A( J. I ) =HOLD 


011500 




INTERCHANGE COLUMNS 


011510 


35 


I *M( X ) 


011520 




IFIM(X)-X) 45.45.37 


011530 


37 


D040 J* 1 »N 


011540 




HOLD=-A( J.K) 


011550 




A( J.X) =A( J, I ) 


011560 


40 


A( J. I ) “HOLD 


011570 




DIVIDE COLUMN BY MINUS PIVOT 


011580 


45 


D055 1=1. N 


011590 


46 


I F ( I-X) 50.55.50 


011600 


50 


A( I .X) =A( I .X ) /( -A(X.X) ) 


011610 


55 


CONTINUE 


011620 




REDUCE MATRIX 


011630 




D065 1*1. N 


011640 




D065 J=1.N 


011650 


56 


IF(I-X) 57,65.57 


011660 


57 


IF(J-IC) 60.65.60 


011670 


60 


AII.J)=A(I»K)*A(K.J)+A(I.J) 


011680 


65 


CONTINUE 


011690 




DIVIDE ROW BY PIVOT 


011700 




D075 J= 1 » N 


011710 


68 


IF( J-IC) 70,75.70 


011720 


70 


A(X»J)=A( K. JI/AIX.K) 


011730 


75 


CONTINUE 


011740 




CONTINUED PRODUCT OF PIVOTS 


011750 



102 



n n 



C 



D = D* A ( K * K ) 

REPLACE PIVOT BY RECIPROCAL 
A(K,K)=1#0/A(K*K) 

80 CONTINUE 
C FINAL ROW AND COLUMN INTERCHANGE 

K = N 

100 K= ( K-l ) 

I F ( K ) 150*150* 103 
103 I = L ( K ) 

I F ( I -K ) 120*120*105 
105 D0110 J = 1 * N 
HOLD = A ( J * K ) 

A ( J *K ) = -A ( J* I ) 

110 A ( J * I ) = HOLD 
120 J=M ( K ) 

IFtJ-K) 100,100*125 
125 D0130 I = 1 * N 
HOLD = A ( K * I ) 

A ( K • I ) =-A ( J# I ) 

130 A ( J , I ) =HOLD 
GO TO 100 
150 RETURN 
END 



C 

C 



SUBROUTINE HANGER ( C LX *CLY * CT7 * U ) 
DIMENSION U ( 6 *6 ) 

DO 1 J S 1 * 6 
DO 1 K= 1 * 6 
1 Ut J*K)=0.0 
UC 1 #1)=1*0 
U(2*2)=1.0 
U(3*3)=l*0 
U( 4,4) -1*0 
U( 5,5) *1 .0 
U( 6 *6 ) =1 *0 
U ( 4 * 3 ) =CTZ 
U( 5,2)=-CLY 
U ( 6 * 1 ) c CLX 
RETURN 
END 



SUBROUTINE HANGEO ( CTX *CLZ *CT Y *U ) 
DIMENSION U t 6 *6 ) 

DO 1 J- 1 # 6 
DO 1 K*l,6 
1 U(J*K)=0.0 
UC 1,1)*1.0 
U ( 2 * 1 ) =CTX 
U( 2 *2) =1*0 
U(3*3)=1.0 
U(4*4)=1.0 
U( 5*4)*CTY 
U ( 5 * 5 ) = 1 *0 
U ( 6 * 3 ) =-CLZ 



011760 

011770 

011780 

011790 

011800 

011.810 

011820 

011830 

011840 

011850 

011860 

011870 

011880 

011890 

011900 

011910 

011920 

011930 

011940 

011950 

011960 

011970 

011980 

011990 

012000 

0*2010 

012020 

012030 

012040 

012050 

012060 

012070 

012080 

012090 

012100 

012110 

012120 

012130 

012140 

012150 

012160 

012170 

012180 

012190 

012200 

012210 

012220 

012230 

012240 

012250 

012260 

012270 

012280 

012290 

012300 

012310 



n n 



U(6.6)=1.0 

RETURN 

END 



SUBROUTINE CFIELO (AJT.R.Z.G.SDRI .RHO.THETA.D.DI • FF AC » A ) 
DIMENSION A ( 6 • 6 ) 

PH I *THETA/Z 
RHOD=RHO 

IF(RHO-0«0)ll»lltl2 

11 V=-1.0 
GO TO 1 

12 V s 1 • 0 

1 CP*COSF ( PH I ) 

SP-S I NF ( PH I ) 

W=1.0/( G * AJT ) 

F1=(PHI*CP + SP) /2 • 0 
F3=(SP-CP*PHI ) /2 • 0 

F5 = ( 2.0-2.0*CP-PHI*SP ) /2.0 
F6 = ( 2.0*PHI+PHI*CP-3.0*SP ) / 2 • 0 
RHO=ABSF(RHO) 

DO 2 J= 1 • 6 
DO 2 K>1.6 

2 A(J.K)=0.0 

I F ( SDR I ) 1 3 . 13 . 14 

13 AREA=( 3.141592654/4.0 ) *(D**2-DI**2 )*FFAC 
SD=RHO*THET A/(G*Z*AREA) 

GO TO 15 

14 SD=0.0 

15 A ( 1 » 1 ) c CP 

A( 1 .2) = (W*RHO*Fl)-(RHO*F3*R) 

A ( 1 • 4 ) = V*SP 

A < 1 . 5 ) = V* ( W+R ) * ( RHO*PH I *SP ) / 2 . 0 
A ( 1 . 6 ) * V* ( W+R ) * ( RHO**2 ) *F3 
A ( 2 • 2 ) =CP 
A ( 2 » 5 ) * V*SP 
A ( 2 • 6 ) * V*RHO* ( 1 .0-CP ) 

A ( 3 * 1 ) *-A (2*6) 

A<3.2)*-1.0*A<1.6) 

A ( 3.301.0 
A ( 3*4) =RHO*$P 

A ( 3 # 5) * ( ( R* ( RHO**2 )*PH1*SP)/2.0)-W*RHO**2*F5 
A ( 3.60 (R*RHO**3*F3)-< W*F6*RHO** 3 ) -SD 
A ( 4 # 1 ) *-V*SP 

A( 4.2 0-VM W+R) *RHO*PHI*$P/2.0 
A ( 4 . 4 ) =CP 

A ( 4 • 5 O ( R*RHO*F 1 ) - ( W*RHO*F 3 ) 

A ( 4 . 6 O A ( 3.5 ) 

A(5.2)«-V*SP 
A< 5.50CP 
A ( 5 . 6 O A ( 3.4) 

A ( 6 * 6 0 1 .0 
RHO=RHOD 
RETURN 
END 



012320 

012330 

012340 

012350 

012360 

012370 

012380 

012390 

012400 

012410 

012420 

012430 

012440 

012450 

012460 

012470 

012480 

012490 

012500 

012510 

012520 

012530 

012540 

012550 

012560 

012570 

012580 

012590 

012600 

012610 

012620 

012630 

012640 

012650 

012660 

012670 

012680 

012690 

012700 

012710 

012720 

012730 

012740 

012750 

012760 

012770 

012780 

012790 

012800 

012810 

012820 

012830 

012840 

012850 

012860 

012870 



* # 
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SUBROUTINE POINO ( AMU * AI X ,AI Y , OMEGA *TL *Z *SDR I *B ) 


012880 

012890 




DIMENSION B ( 6 *6 ) 


012900 




AM=AMU*TL/( Z-liO) 


012910 




DO 1 J= 1 * 6 


012920 




DO 1 K = 1 * 6 


012930 


1 


B( J*K) =0.0 


012940 




B(l»l)=1.0 


012950 




B ( 2 * 1 ) = -AM* ( A I X*OMEGA ) **2 


012960 




B<2*2)=1.0 


012970 




B<3*3)=1.0 


012980 




B(4*4)=1.0 


012990 




I F ( SDR I -0* ) 9*9. 10 


013000 


9 


B ( 5 * 4 ) =-AM* ( A I Y*OMEGA ) ** 2 


013010 




GO TO 11 


013020 


0 


B ( 5 * 4 ) =0 


013030 


.1 


B<5*5)=1.0 


013040 




B ( 6 * 3 ) = AM*OMEGA**2 


013050 




B(6*6)=1.0 


013060 




RETURN 


013070 




END 


013080 






013090 

013100 




SUBROUTINE CF I ELD (R*Z*G*SD* RHO *THETA*D*DI*E*A) 


013110 




DIMENSION A ( 6 * 6 ) 


013120 




PH I =THETA/Z 


013130 




RD= ABSF ( RHO ) 


013140 




I F ( RHO ) 1 * 1 * 2 


013150 


1 


V=-1.0 


013160 




GO TO 3 


013170 


2 


V= 1 • 0 


013180 


3 


CP=COSF(PHI ) 


013190 




SP=S I NF (PHI ) 


013200 




DO 4 J= 1 * 6 


013210 




DO 4 K= 1 * 6 


013220 


4 


A< J*K)=0.0 


013230 




A ( 1*1)=CP 


013240 




A ( 1*2) --V*SP 


013250 




A ( 1 ,3)=-V*RD*( 1.0-CP) 


013260 




A ( 1 ,4) =-V*RD**2*R* (PHI-SP) 


013270 




AR=3.141592654*(D**2-DI**2 )/4.0 • 


013280 




W = 1 • 0/ ( E* AR ) 


013290 




F3= ( SP-PH I *CP ) * • 5 


013300 




F1=(PHI*CP+SP)*.5 


013310 




F5= <2.0-2*0*CP-PHI*SP) *.5 


013320 




F6= (2.0*PHI+PHI*CP-3.0*SP )**5 


013330 




A ( 1 .5) =V*(RD*W*PHI*SP*.5-R*F5*RD**3) 


013340 




A( 1 *6)=RD*W*F1+R*F6*RD**3 


013350 




A ( 2 • 1 ) =V*SP 


013360 




A ( 2 * 2 ) =CP 


013370 




A ( 2 * 3 ) *RD* SP 


013380 




A( 2*4) =R*RD**2* (1.0-CP ) 


013390 




A( 2*5) =RD*F3*( W+R*RD**2 ) 


013400 




A ( 2 * 6 ) = A ( 1*5) 


013410 




A(3*3)=1.0 


013420 




A( 3*4)=RD*R*PHI 


013430 



nr>nnn n n 



A< 3»5)=A(2.4) 
A ( 3 • 6 ) = A ( 1*4) 
A(4»4)=1.0 
A(4,5)=A(2.3) 
A ( 4 » 6 ) = A ( 1*3) 
A ( 5 • 5 ) = CP 
A(5.6)=A( U2 ) 
A(6*5)=A(2.1 ) 
A ( 6 * 6 ) 3 CP 
RETURN 
END 



SUBROUTINE POINT < AMU . A I Y * OMEGA * TL • Z 
DIMENSION B ( 6 »6 ) 

AM = AMU* TL / ( 2- 1 • ) 

DO 1 J= 1 * 6 
DO 1 K = 1 * 6 

1 B ( J * K ) =0 • 0 
B(ltl)-1#0 
B(2»2)*l*0 
B(3t3)=1.0 

I F ( R I ) 2 * 2 1 3 

2 B(4t3)*-AM*(AIY*OMEGA)**2 

3 B(4 t 4)*1.0 

B( 5 *2 ) =AM*OMEGA**2 

B<5t5)=1.0 

B ( 6 • 1 ) =-B ( 5 * 2 ) 

B ( 6 1 6 ) * 1 • 0 

RETURN 

END 



SUBROUTINE SF I ELD < D I • TL * R , 2 ,G . SD • D • 
DIMENSION A ( 6 • 6 ) 

ARA=( 3# 141592654/4.0) *<D**2-Dl**2) 
AREA=ARA*FFAC 
DL =TL/2 
DO 1 J*2 * 5 

1 A(ltJ)*Q.O 
A(ltl)*1.0 

A( 1 • 6 ) *DL/ ( E*ARA ) 

A { 2 * 1 ) *0. 0 
A(2#2)«1.0 
A ( 2 • 3 ) *DL 

A ( 2 .4) *DL**2*R/2.0 
I F ( SD ) 2 *2 • 3 

2 A(2t5)=DL**3*R/6.0-DL/ (G*AREA ) 

GO TO 4 

3 A(2t5)»DL**3*R/6.0 

4 A( 2 • 6 ) =0.0 
A ( 3 * 1 ) *0.0 
A ( 3 t 2 ) *0 . 0 



013440 

013450 

013460 

013470 

013480 

013490 

013500 

013510 

013520 

013530 

013540 

013550 

013560 

RI.B) 013570 

013580 

013590 

013600 

013610 

013620 

013630 

013640 

013650 

013660 

013670 

013680 

013690 

013700 

013710 

013720 

013730 

013740 

013750 

013760 

013770 

013780 

013790 

’ • F F A C • A ) 013800 

013810 

013820 

013830 

013840 

013850 

013860 

013870 

013880 

013890 

013900 

013910 

013920 

013930 

013940 

013950 

013960 

013970 

013980 

013990 
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01A000 


A(3*3>=1.0 


01A010 


A ( 3 * A ) = DL*R 


01A020 


A<3»5)=A(2*A) 


01A030 


A ( 3.6) *0.0 


01A0A0 


DO 5 J = 1 • 6 


01A050 


A ( A * J ) =0.0 


01A060 


A ( 5 * J ) =0.0 


01A070 


, A ( 6 * J ) =0.0 


01A080 


A( A, A) =1.0 


01A090 


A ( A . 5 ) =DL 


01 A 100 


A(5. 51=1.0 


01A110 


A( 6 *6 ) = 1 .0 


01A120 


RETURN 


01A130 


END 


01A1A0 

01A150 


SUBROUTINE SFIELO ( D I .TL . A JT .R . Z .0 .SDR I . D. FF AC . A ) 


01A160 

01A170 


DIMENSION A ( 6 . 6 ) 


01A180 


AREA=<3.1A159265A/A.0)*(D**2-DI*»2)*FFAC 


01A190 


dl=tl/z 


01A200 


DO 2 J = 1 . 6 


01A210 


DO 2 <=1.6 


01A220 


2 A ( J .< ) “0. 0 


01A230 


All. 11=1.0 


01A2A0 


A(2.2)*1.0 


01A250 


A(3.3)=1.0 


01A260 


A ( A .A) *1 .0 


01A270 


A(5.5)*1.0 


01A280 


A(6.6)=1.0 


01A290 


IFISDRI-O. 17.7.8 


01A300 


7 A( 3.6) =DL*+3*R/6.0-DL/ (G*AREA) 


01A310 


GO TO 1 


01A320 


8 A(3.6)=DL»*3*R/6.0 


01A330 


1 A(1.2)=DL/(G*AJT) 


01A3A0 


A ( 3 .A ) =DL 


01A350 


A(3.5)=DL**2*R/2.0 


01A360 


A ( A .5) =DL*R 


01A370 


A(A.6)=A(3.5) 


01A380 


A( 5 .6 ) =DL 


01 A 390 


RETURN 


01AA00 


END 


01AA10 

01AA20 


SUBROUTINE MATMUL (A.B.U) 


01AA30 

01AAA0 


DIMENSION A(6.6).B(6.6).C(6.6) ,U(6.6) 


01AA50 


DO 1 J=1 »6 


01AA60 


DO 1 <=1*6 


01AA70 


C( J. 0=0.0 


01AA80 


DO 1 L=1 *6 


01AA90 


1 C( J.O=C( J.O+B(J.L)*A(L.<) 


01A500 


DO 2 J = 1 .6 


01A510 


DO 2 <=1.6 


01A520 


2 U< J.O=C< J.O 


01A530 


RETURN 


01A5A0 


END 


01A550 



\o r i 



n n 





SUBROUTINE F I NBR A ( U * VV * V * MAT * A * Z ) 


014560 

014570 




DIMENSION U C 6 * 6 ) * A C 6 • 6 ) *VV ( 6 * 3 ) *V( 6 • 3 ) 


014580 




TYPE DOUBLE VV*V 


014590 




I F ( MAT ) 65*64*65 


014600 


64 


DO 67 J = 1 * 6 


014*10 




DO 67 K = 1 • 3 


014620 




V( J*K)=0. 


014630 




DO 67 L = l*6 


014640 


67 


V( J.K)=V( J.K)+UCJ*L)*VV(L*K) 


014650 




RETURN 


014660 


65 


N=2-l. 


014670 




DO 5 M»1*N 


014680 




DO 4 J= 1 * 6 


014690 




DO 4 K= 1 • 3 


014700 




V( J.O *0. 


014710 




DO 4 L * 1 • 6 


014720 


4 


V< J*K)*V( J * K ) -KU ( J*L)*VV(L*K) 


014730 




DO 5 J r 1 * 6 


014740 




DO 5 K s 1 * 3 


014750 


5 


VV( J*K)*V( J*K) 


014760 




DO 7 J = 1 * 6 


014770 




DO 7 K = 1 *3 


014780 




V< J*K)=0* 


014790 




DO 7 L * 1 * 6 


014800 


7 


V(J»K)sV(J*K)+A(J*L)*VV(L»K) 


014810 




RETURN 


014020 




END 


014030 






014840 

014850 




SUBROUTINE FINMAT ( U * UU * V . MAT * A . Z ) 


014860 




DIMENSION U ( 6 *6 ) * A ( 6 . 6 ) * V ( 6 • 3 ) *UU ( 6 * 3 ) 


014870 




TYPE DOUBLE UU*V 


014880 




I F ( MAT ) 65*64*65 


014890 


64 


DO 67 J~1 *6 


014900 




DO 67 K-l * 3 


014910 




V( J*K)*0. 


014920 




DO 67 L = 1 * 6 


014930 


67 


V( J*K)*V< J.K)+U( J»L)*UU(L*K) 


014940 




RETURN 


014950 


65 


N=Z- 1 • 


014960 




DO 5 M= 1 * N 


014970 




DO 4 J = 1 * 6 


014980 




DO 4 K 3 1 * 3 


014990 




V< J*K)*0# 


015000 




DO 4 L-l *6 


015010 


4 


V(J*K)=V(J*K)+UIJ*L) *UU C L * K ) 


015020 




DO 5 J -1 * 6 


015030 




DO 5 K = 1 ♦ 3 


015040 


5 


UU(J*K)*V(J..K) 


015050 




DO 7 J« 1 * 6 


015060 




DO 7 K * 1 • 3 


015070 




V( J*K)=0. 


015080 




DO 7 L - 1 * 6 


015090 


7 


V(J*K) 3 V(J*K)+A(J*L)*UU(L*K) 


015100 




RETURN 


015110 



I 08 



c 

c 



END 



SUBROUTINE BRANCH ( R 3 » VV ♦ PH I * U ) 

DIMENSION R 3 ( 3 • 3 ) ♦ VV (6*3) * U ( 6 * 6 ) *Rl(25t25) * L L ( 2 5 ) *MM (25)# R2 (3*3)# 
1 R ( 3 ♦ 3 ) f G1 ( 3 « 3 ) «G2 ( 3 * 3 ) 

TYPE DOUBLE VV 
DO A L=1.3 
R1(1*L)=VV( ltL) • 

R1(2'L)«VV(2*L) 

R 1 ( 3 ♦ L ) =VV ( 3 * L ) 

R2( 1 *L ) =VV ( 4 ♦ L ) 

R2(2*L)=VV(5'*L) 

A R2I 3 ♦ L ) =VV ( 6 * L ) 

SP = S I NF ( PHI ) 

CP=COSF(PHI ) 

G 1 ( 1 # 1 ) =CP 
G 1 ( 1 • 2 ) *-SP 
G1 ( 1 1 3 ) =0 • 

G 1 ( 2 • 1 ) =SP 
G 1 ( 2 1 2 ) =CP 
G 1 ( 2 » 3 ) = 0 • 

G 1 C 3 ♦ 1 ) *0 • 

Gl(3#2)=0. 

G1 ( 3 1 3 ) *1 • 

G2 C 1 ♦ 1 )=1. 

G2 ( 1 »2 ) =0. 

G2 ( 1 • 3 ) =0 • 

G2 ( 2 1 1 ) =0 • 

G2 ( 2 * 2 ) =CP 
G2 ( 2 * 3 ) s -SP 
G2 ( 3 ♦ 1 ) =0 • 

G2 ( 3 f 2 ) =SP 
G2 ( 3 • 3 ) =CP 

CALL INVERT ( R 1 ♦ 3 ♦ Dt LL *MM ) 

DO 5 J = 1 *3 
DO 5 K = 1 « 3 
R3( J*K)=0.0 
DO 5 L= 1 • 3 

5R3(JfK)=R3(J*K)+Rl(J»L)*Gl(L*K) 

DO 6 J= 1 • 3 
DO 6 K= 1 » 3 
Rl( J*K)=0.0 
DO 6 L = 1 * 3 

6 R1 C J*K) *R1( J«K)+G2C JtL)»R2(L*K) 

DO 7 J= 1 ♦ 3 

DO 7 K= 1 ♦ 3 
R( JfK)=0.0 
DO 7 L = l*3 

7 R(J*K)=R(J#K')+R1(J#L)*R31L*K) 

DO 8 J*l#6 

DO 8 K*l*6 

8 U(J*K)=0.0 
U( 1*1)*1.0 
U ( 2 ♦ 2 ) = 1 *0 
U ( 3 • 3 ) =1 #0 



015120 
015130 
0 15 1 AO 
01 51 50 
015160 
015170 
015180 
015190 
015200 
015210 
015220 
015230 
015240 
015250 
015260 
015270 
015280 
015290 
015300 
015310 
015320 
015330 
015340 
015350 
015360 
015370 
015380 
015390 
015400 
015410 
015420 
015430 
015440 
015450 
015460 
015470 
015480 
015490 
015500 
015510 
015520 
015530 
015540 
015550 
015560 
015570 
015580 
015590 
015600 
015610 
015620 
015630 
015640 
015650 
015660 
015670 



I 0? 



U( 4*4) =1.0 


015680 


U( 5*5) =1 .0 


015690 


U( 6*6 ) 3 1 *0 


015700 


DO 9 L 3 1 » 3 


015710 


U(4»L)=R( 1*L) 


015720 


U(5*L)=R(2*L) 


015730 


U(6*L)=R(3*L) 


015740 


RETURN 


015750 


END 


015760 




015770 




015700 


SUBROUTINE BRANCO ( R 3 * VV* PH I *U ) 


015790 


DIMENSION R3( 3*3) *U(6*6) * VV ( 6 * 3 ) *R1 (25*25) »R2(3*3) • 


015800 


TYPE DOUBLE VV 


015010 


1G1(3*3) *G2(3*3) *S(3*3) *LL(25)* MM (25) 


015020 


DO 4 K* 1 * 3 


015030 


Rl(ltK) 3 VV( 1 * K ) 


015040 


R 1 ( 2 »K ) *VV ( 3 * K ) 


015050 


R1(3*K)*VV(4*K) 


015060 


R2 ( 1 »K ) 3 VV ( 2 * K ) 


015070 


R2 ( 2 * K ) =VV ( 5 * K ) 


015000 


R2 ( 3 *K ) 3 VV ( 6 * K ) 


015890 


SP=SINF(PHI ) 


015900 


CP=COSF ( PH I ) 


015910 


G1 ( 1 * 1 ) =CP 


015920 


G1 ( 1 * 2 ) =0 • 


015930 


G1(1»3)«-SP 


015940 


G 1 ( 2 * 1 ) =0 • 


015950 


G1 ( 2 *2 ) =1 • 


015960 


G 1 ( 2 1 3 ) =0 • 


015970 


G1 (3*1) *SP 


015980 


G1 ( 3 * 2 ) =0 • 


015990 


G1 ( 3 * 3 ) 3 CP 


016000 


G2 ( 1 * 1 ) 3 CP 


016010 


G2 ( 1 * 2 ) 3 SP 


016020 


G2 ( 1 * 3 ) *0 • 


016030 


G2 ( 2 * 1 ) =-SP 


016040 


G2 ( 2 *2 ) 3 CP 


016050 


G2 ( 2 * 3 ) 3 0 • 


016060 


G2 ( 3 • 1 ) =0 • 


016070 


G2 ( 3 * 2 ) =0 • 


016000 


G2 ( 3 * 3 ) = 1 • 


016090 


CALL INVERT ( R1 * 3 * D. LL *MM ) 


016100 


DO 5 J = 1 * 3 


016110 


DO 5 K = 1 * 3 


016120 


R3( J*K)=0.0 


016130 


DO 5 L= 1 ♦ 3 


016140 


> R3( J*K)=R3( J*K)+R1< J»L)*G1 (L*K) 


016150 


DO 6 J 3 1 * 3 


016160 


DO .6 K 3 1 * 3 


016170 


R 1 ( J*K)=0*0 


016100 


DO 6 L 3 1 * 3 


016190 


► R1(J*K)=R1(J*K)+G2(J*L)*R2(L*K) 


016200 


DO 7 J= 1 * 3 


016210 


DO 7 K= 1 * 3 


016220 


S( J*M=0.0 


016230 



! I 0 
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DO 7 L= 1 ♦ 3 

7 S(MI=5IJ.K| t RllJ.L)*R3(L,K) 

DO 8 J = l*6 
DO 8 <*1*6 

8 U(J*<)=0.0 

U(2.1)=S(1.1) 

U«2*2)=1.0 

U(2*3) =5(1*2) 

U( 2 *4 ) =S I 1*3) 

U<3.3)=1.0 

U<4*4)=1.0 

U(5*1)=S(2*1) 

U(5*3)=S(2*2) 

U(5*4)=S(2«3) 

U< 5 * 5 ) = 1 .0 
U( 6 * 1 ) =S( 3 . 1 ) 

U(6.3)»S(3.2) 

U(6*4)=S(3*3) 

U(6.6) =1.0 

RETURN 

END 

SUBROUTINE ST ATEM ( SV I *UU* D1 *02 . PI * P2 • PP * PM * ss > 
DIMENSION SVI 16*1 ) *UU(6*3) 

TYPE DOUBLE UU« PI * P2 *D1 • 02 
M = 0 

DO 1 J = 1*6 
DO 1 <*1*3 
1 UU(J*<)=0. 



016240 

016250 

016260 

016270 

016280 

016290 

016300 

016310 

016320 

016330 

016340 

016350 

016360 

016370 

016380 

016390 

016400 

016410 

016420 

016430 

016440 

016450 

016460 

016470 

016480 

016490 

016500 

016510 

016520 

016530 

016540 

016550 



IF(PM) 40*40*41 


016560 


40 DO 42 J = l*6 


016570 


I F < 5V I ( J * 1 ) ) 42*42*43 


016580 


43 M=M+ 1 


016590 


UU( J*M)=1« 


016600 


IF(M-l) 42*44*42 


016610 


44 UU ( J * 2 ) =-D 1 


016620 


UU< J*3 ) =-02 


016630 


42 CONTINUE 


016640 


RETURN 


016650 


41 PI =P 1+0 1 


016660 


P2=P2+02 


016670 


IF(SS) 14*15.15 


016680 


15 IF ( SS- 1 • ) 12*11*10 


016690 


14 I F ( PP- 1 • ) 12*11*10 


016700 


10 DO 4 J = 1 • 6 


016710 


I F ( SV I < J*l) > 4*4*2 


016720 


2 M=M+ 1 


016730 


UU(J*M)=1. 


016740 


I F ( M-2 ) 4.3.5 


016750 


3 UU(J*1)-P1 


016760 


GO TO 4 


016770 


5 UU ( J • 1 ) =P2 


016780 


4 CONTINUE 





I I I 



1 1 .DO 21 J = 1 *6 

IF ( SVI ( J t 1 ) ) 21*21*22 

2 2 , / = i v i+ 1 

UU( J* v «) =1. 

IF(M-l) 25*26*25 
26 UU< J*3) =P1 
2 6 IF(M-2) 21 *26*21 

2 3 him J *2 ) =D? 

21 rONTINUF 
P FT URN 

12 DO 21 J = 1 * 6 

IF* SVI ( J * 1 ) ) 31*31*32 
32 V« = M+1 

UU ( J *M ) = 1 • 

IF(M-l) 35*36*35 
36 UU ( J * 2 ) =P2 
35 [F(M-2) 31,^3*33 

3 3 UU ( J * 2 ) =P 1 
3 1 CONTIN'JF 

RFfuPN 

CMH 

r 

r 

9UBR0UT INF STATFB* SV** VV*N ,YY,L*BRC 
DIMENSION VV(6*3) »YY( 3,3*25) • V ( 3 * 3 
TYPF DOUBLE VV 
31=0. 

B2 = 0 • 

B3 = 0 • 

D 1 = 0 « 

D2 = 0. 

IF (PM) 33,40,22 

22 IF<PRM 21*20*21 

23 I f ( L- 1 ) 23,20,23 

23 DO )0‘ J=1 ,3 

DO 10 K = 1 

1C V ( J * K ) = Y Y ( J , < , N ) 

IF(SS) 11*12*12 

11 I F ( PP- 1 • ) 32*31,30 

12 IF<SS-1.) 33*31,30 

3 ' 1 Bl=V(l*l)+V( 1*2)*D1+V(1*3) *D2 + D l 
92 = V ( 2 * 1 ) +V ( 2 * 2 ) *f) 1 +V ( 2 » 3 ) *D2+B2 
rt = V( 3 , 1 ) +v ( 3 ,2 ) *D1+V< ° , 3 ) *02 + 03 
00 TO 33 

?. 1 R1 =«1 + V( 1 * 1 )*D 2 +y ( 1 , ) #pi + v( 1 , ? ) 

P2 = F2+V( 2 , 1 ) *D?+V ( 3 , 3 ) *D1 +V( 2 , ? ) 
i^ = R3+V( 3 , 11 *D 2+V( * , 3 ) *U1+V< ? , 2 ) 

00 TO 33 

32 Bl=31+V< 1 * 1 ) *D1+V( 1 * 2 ) 02 + V ( 1 , 3 ) 

F2 = B2+V(2 * 1 ) Ol+V ( 2,2 ) *D2+V( 2 ) 

P3 = *3 + V( 3* 1 )*D1+V( 3*2 ) *D2 + V( 3,3 ) 

3 3 01= 0 1 *■ R 1 

02=Bl*02+c2 
c* = « 1 -^ 03-f.n -a 



01679Q 

016800 

016810 

016820 

016830 

016840 

016850 

016860 

016870 

016880 

016890 

016900 

016910 

016920 

016930 

016940 

016950 

016960 

016970 

016980 

016990 

017OO0 

017010 

017020 

01 ,02, PP, PM, S3, Q1,Q2*03, F INK) 017030 

,SV6(6*22) 017040 

017050 

017060 

017070 

017080 

017090 

017100 

017110 

01712" 

017130 

017140 

017150 

017160 

017170 

017180 

017190 

017200 

017210 

017220 

017230 

017240 

017250 

017260 

017270 

017280 

017290 

017300 

017310 

017320 

017330 



r> r> 



20 


M= 1 


017340 




DO 1 J = 1 ,6 


017350 




DO 1 K = 1 , 3 


017360 


1 


VV( J, 0=0.0 


017370 




DO 4 J= 1 » 6 


017380 




I F ( S\/B < J ,N ) ) 4,4,2 


017390 


2 


I F ( M— 1 ) 7 , 3 , 7 


017400 


3 


I F ( PM ) 34,35,35 


017410 


35 


VV< Jil)*01 


017420 




GO TO 36 


017430 


34 


VV(J,1)=1. 


017440 


36 


M=M+1 


017450 




GO TO 4 


017460 


7 


I F { M-2 ) 4,8,9 


017470 


8 


W(J,1) =Q2 


017480 




W ( J »M ) *1 • 


017490 




M = M+ 1 


017500 




GO TO 4 


017510 


9 


VV< J»1)=Q3 


017520 




VVC J,M)=1. 


017530 


4 


CONTINUE 


017540 




RETURN 


017550 


40 


IFCFINK) 54,51,54 


017560 


54 


IF(BRC) 51,51,52 


017570 


52 


DO 53 J = 1 , 3 


017580 




DO 53 < = 1,3 


017590 


5 3 


V< J,<)*YY< J»K»N) 


017600 




D1=(V(1»2)/V( 1 , 1 ) +V( 2 ,2 ) / V( 2 • 1 ) +V< 3,2 ) / V< 3, 1 ) ) / 3 • 


017610 




D2=IV(1»3)/V(1»1)+V<2»3)/\M2,1)+V(3,3)/V<3,1 ) )/3# 


017620 


51 


M = 0 


017630 




DO 41 J = l,6 


017640 




DO 41 <=1,3 


017650 


41 


W< J»<)*0. 


017660 




DO 42 J = 1 ,6 


017670 




I F { SVB ( J , N ) ) 42,42,43 


017680 


43 


M*M+1 


017690 




VV< J,M) =1. 


017700 




IF<M-1) 42,44,4 2 


017710 


44 


VV( J,2)*-D1 


017720 




VV( J,3)=-D2 


017730 


42 


CONTINUE 


017740 




RETURN 


017750 




END 


017760 

017770 

017780 




SUBROUTINE BRCORI R3,UU»XX* J<, VV ,PM ) 


017790 




DIMENSION R 3 1 3,3) , UU ( 6 ,3 ) , XX ( 3 » 3 ) ,UUU< 3*3) ,VV<6»3 ) 


017800 




TYPE DOUBLE VV,UU 


017810 




IF(PM) 20*21,20 


017820 


20 


I F < J< ) 10,10,11 


017830 


10 


DO 12 J=1 , 3 


017840 




DO 12 <=1,3 


017850 


12 


UUU( J,<)=UUU,<) 


017860 




GO TO 3 


017870 


11 


DO 13 <=1,3 


017880 




UUU( 1,K)=UU( 1,0 


017890 



I 13 



n n 





UUU<2,K)*UU< 3.10 


017900 


13 


UUU< 3,K.)=UU<4.K> 


017910 


3 


00 1 J* 1 » 3 


017920 




00 1 K = 1 * 3 


017930 




XX( 


017940 




00 1 L« 1 * 3 


017950 


1 


XX ( J.K)»XX( J,K)+R3( J.L >*UUU(L,IO 


017960 




RETURN 


017970 


21 


IF'(JK) 25,25.22 


017980 


25 


00 23 J = 1 ,3 


017990 




DO 23 K-1.3 


018000 


23 


XX( J.K)«VV( J.K) 


018010 




RETURN 


018020 


22 


00 24 ,3 


018030 




XX( l.K)«VV( 1 > K ) 


018040 




XX ( 2 ,K ) “VV ( 3 »K ) 


018050 


24 


XX ( 3 »K. ) *VV ( 4 » K ) 


018060 




RETURN 


018070 




END 


018080 

018090 

018100 




SUBROUTINE OELMAI UU, SVE , V . X2 , Y 2 .KK.K , X 1 ,FF ,PP . HH , RR , SS » Y 1 , REM ,F I NK. , 


018110 


1PM * D1 * 02 ) 


018120 




DIMENSION UU (6*3) *V(6*3)*SVEI6*1) *0(3*3) 


018130 




TYPE DOUBLE UU * V * 0 1 X 1 , X2 * Y 1 * Y2 , DV * 01 * 02 


018140 




L = 0 


018150 




X1*0. 


018160 




X2«0. 


018170 




Y2*0. 


018180 




Y1*0. 


018190 




D1 *0 • 


018200 




02 = 0* 


018210 




00 111 J* 1 * 3 


018220 




00 111 K* 1 * 3 


018230 


1U 


D( J,K> *0* 


018240 




00 91 N«l*6 


018250 




I F ( SVE ( N * 1 ) ) 91*2*91 


018260 


2 


L*L+ 1 


018270 




DO 92 K = 1 * 3 


018280 


92 


V ( L • K ) *UU ( N * K ) 


018290 


91 


CONTINUE 


018300 




IF(PM) 66*64*61 


018310 


64 


IF(FINK) 63*62*63 


018320 


62 


Dl® ( V( 1 • 2 ) / V C 1 * 1 ) +V( 2*2)/V(2*l)+V( 3*2)/V(3*l ) )/3# 


018330 




D2*(V(l*3)/V(l*l)+V(2*3)/V(2*l)+V(3*3)/V(3*l))/3» 


018340 




FINK«1. 


018350 




RETURN 


018360 


63 


F I NK *0 • 


018370 




RETURN 


018380 


66 


RETURN 


018390 


61 


IF(SS) 31 1 » 31 2 • 31 2 


018400 


311 


IF ( PP— 1 • ) 400 * 40 1*402 


018410 


312 


IF(SS-1«) 313t314 t 315 


018420 


402 


I F ( RR- 1 • ) 300*403*300 


018430 


315 


I F ( RR- 1 • ) 340*341*340 


018440 


300 


I F ( V ( 3 * 3 ) ) 93,94,93 


018450 



I l 4- 



340 I F ( V ( 3 * 2 ) ) 93*94*93 

93 DO 95 J = 2 * 3 
DO 95 K= 1 * 3 

95 D ( J * K ) =* V ( J * K ) 

DV=V(1 *2)*D(2*3)-V( 1 • 3 ) *D < 2 * 2 ) 

IF(DV) 98*94*98 

94 RR= 1*0 
HH= 1 • 

RETURN 

403 I F ( V ( 2 * 3 ) ) 96*200*96 

341 I F I V ( 2 • 2 ) )96*200*96 

96 DO 97 K = 1 * 3 
D(2*K)=V(3*K) 

97 D ( 3 * K ) s V { 2 * K ) 

DV = V( 1.2)*D(2.3)-V( 1 * 3 ) *D ( 2 * 2 ) 

IF(DV) 98*200.98 

98 X1=X1-(V< 1.1)*D(2.3)-V<1.3)*D(2.1> )/DV 
X2=X2- ( V ( 1.2 )*D(2.1)-V(1*1)*D{ 2*2) )/DV 
IF (XI) 100*102*100 

102 I F ( X 2 ) 100.101*100 
101 FF= 1 • 0 
RETURN 

100 IF(SS) 354.355*355 

354 Y2=Y2-D(3*1 )/D( 3.3 )-Xl*D( 3*2 )/D( 3*3) 
RETURN 

355 Y1=Y1— D(3*l )/D( 3.2 ) “X2*D( 3 * 3 ) /D ( 3*2) 

. RETURN 

99 IF(SS) 20*21.21 

20 PP=1. 

GO TO 23 

21 5S= 1 • 

23 RR=0 • 

HH* 1 • 

RETURN 

314 I F ( RR- 1 • ) 330*331.330 
401 IF(RR-l.O) 301*404.301 

330 IF ( V( 3 * 1 ) ) 81*82.81 
301 I F ( V ( 3 * 2 ) ) 81.82.81 

81 DO 83 J = 2 * 3 
DO 83 K = 1 *3 

83 D(J*K)=V(J*K) 

DV=V< 1 .1 ) *D ( 2 .2 )-V< 1 *2 )*D( 2.1 ) 

IF(DV) 86.82*86 

82 RR=1.0 
HH * 1 • 

RETURN 

331 I F ( V ( 2 * 1 ) ) 84.70*84 

404 I F ( V ( 2 • 2 ) ) 84.70*84 

84 DO 85 K=1 *3 
D(2*K)«V(3»K) 

85 D(3*K)=V(2*K) 
DV*VU.1)*D(2.2)-V(1.2)*D(2*1) 

IF(DV) 86*70.86 

8 6 X1=X1-(V( 1*3)*D(2*2) -V <1.2 )*D( 2*3) )/DV 
X2-X2-(V( 1*1 ) *D ( 2 * 3 ) -V (1*3)*D(2*1) )/OV 
IF(X1)2C3. 201*203 



018460 

018470 

018480 

018490 

018500 

018510 

018520 

018530 

018540 

018550 

018560 

018570 

018580 

018590 

018600 

018610 

018620 

018630 

018640 

018650 

018660 

018670 

018680 

018690 

018700 

018710 

018720 

018730 

018740 

018750 

018760 

018770 

018780 

018790 

018800 

018810 

018820 

018830 

018840 

018850 

018860 

018870 

018880 

018890 

018900 

018910 

018920 

018930 

018940 

018950 

018960 

018970 

018980 

018990 

019000 

019010 



?% 



201 


IF ( X 2 > 203*202*203 


019020 


202 


FF=1.0 


019030 




RETURN 


019040 


203 


IF(SS) 352.353*353 


019050 


353 


Yl=Yl-D(3 .3 ) /D( 3 * 1 ) -X2*D< 3 *2 ) /D< 3 * 1) 


019060 




RETURN 


019070 


352 


Y2 = Y2-D(3*3 )/D( 3.2 ) -X 1*D( 3 . 1 >7D< 3 . 2 ) 


019Q80 




RETURN 


019090 


70 


IF(SS) 25*24*24 


019100 


24 


SS = 2 • 


019110 




GO TO 26 


019120 


25 


PP = 2. 


019130 


26 


RR = 0 • 


019140 




HH = 1 • 


019150 




RETURN 


019160 


313 


IFIRR-1. ) 320.321.320 


019170 


400 


IFIRR-1.) 302*405.302 


019180 


320 


IF CV(3*3) ) 71.72.71 


019190 


302 


IF ( V ( 3 . 1 ) ) 7 1 * 72 • 7 1 


019200 


71 


DO 73 J-2 • 3 


019210 




DO 73 K= 1 * 3 


019220 


73 


D( J.K)=V( J.K) 


019230 




DV= V (1*1) *D( 2 .3 )-V( 1 *3 ) *D( 2 .1 ) 


019240 




IF(DV) 76.72.76 


019250 


72 


RR=1.0 


019260 




HH = 1 • 


019270 




RETURN 


019280 


40 5‘ 


I F ( V ( 2 * 1 ) ) 74.99.74 


019290 


321 


I F ( V ( 2 . 3 ) ) 74.99.74 


019300 


200 


KKK*1 


019310 




RR = 0 • 


019320 




RETURN 


019330 


74 


DO 75 K = 1 . 3 


019340 




D(2.K)*V(3.K) 


019350 


75 


D(3.K)«V(2.K) 


019360 




DV*V ( 1 . 1 ) *D C 2 . 3 ) -V ( 1 . 3 ) *D ( 2 . 1 ) 


019370 




IF(DV) 76.99.76 


019380 


76 


X2 -X2- ( V C 1 . 2 ) # D(2*3)-V(1*3) # D(2*2) )/OV 


019390 




X1«X1-(V( 1.1 )*D(2.2)-D(2.1 )*V(1 .2) )/DV 


019400 




I F ( X 1 ) 213.211.213 


019410 


211 


IF ( X 2 ) 213.212.213 


019420 


212 


FF = 1 • 0 


019430 




RETURN 


019440 


213 


IF(SS) 350.351.351 


019450 


351 


Y1«Y1-CD(3.2)+X2*D«3.1 ) )/D(3*3) 


019460 




RETURN 


019470 


350 


Y2 e Y2-(D(3.2)+Xl # D(3.3) )/D(3.1) 


019480 




RETURN 


019490 




END 


019500 



APPENDIX E 



ACCURACY OF METHODS & SAMPLE PROBLEMS 
E-l General remarks. 

Appendix E contains the results of vibration analyses of several 
piping systems. The first group was performed to establish the accuracy 
of the methods and compare them with the VIPIPE output. The second group 
was performed to establish assurance of solution by comparing the natural 
frequencies got from various methods. 

E-2 Method accuracy analyses. 

E-2-1 Straight sections. 

End conditions. 

System 1 - fixed-fixed. 

System 2 - fixed-free. 

System 3 - fixed-propped. 

System parameters (for system 1,2,3). 
length - 200 inches, 
diameter - 2.0 inches, 
wall thickness - 0.125 inches. 

Intensive properties (for system 1,2,3) 
density - 490 lbs/cu-ft 
shear modulus - 12x10 psi 
elastic modulus - 30x10 psi 

Mathematical model: Distributed mass with the effect of shear 

deflection and rotational inertia neglected to permit comparison with 
classical solution. 



117 



Comparison of results. 

As explained in Appendix B-9-4, when the iteration reaches the 
solution acceptability criterion, the corresponding natural frequency 
is determined using linear interpolation. Theoretically exact compari- 
son values are obtained by use of the same iteration method with the 
solution acceptability criterion as 0.0007 rad/sec. The values obtained 
using VIPIPE have a solution acceptability criterion 0.0185 rad/sec. 
Comparison is made by percent difference. 

The M-D method fails the systems 1 through 3, because the purging 
factors turn out to be all zeros. So for the purpose of testing the 
M-D method, modified System 1 is provided below. 

The letters D or S in parentheses denote double precision and single 
pteoision respectively. 



Modified System 1: 

Add a very small curved section to the System 1. Second section is 
a mathematical model which has the same property as System 1 except the 
system parameters. 



Second section parameters: 

included angle of arc - 45 degrees 
radius of curvature - 0.01 inches 
(This value of radius of curvature is 
not possible in an actual system, since the radius of curvature less than 
the radius of the pipe. In effect subroutine STIFCO or STIFOO will pro- 
vide the transfer matrix of a massless rigid corner). 
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Results and comparison.* 



a. System 1 and modified System 1.** (in-plane). 





Mode 




Freq. (rad/sec) and difference in %. 


i 




No. 


Comparison values 


Delta m. (D) 


P-M.-A(D) 


P-Mm.-F(D) 




1 


75.09994540 


75.09994544 
0 7. 


75.09994550 
0 % 


r 

75.09994571 
0 %| 




2 


207.01589137 


207.01589036 

5.0xl0' ? 7. 


207.01589094 
0 7 , 


207.01589228 
4.4xl0" 7 7. 




3 


405.83391916 


405.83393183 
3. 0xl0 _6 % 


405.83403022 

2.4x10 -5 % 


405.83394702 
7 .0xl0~ 7 % 




4 


670.86408383 


670.86310856 
2. 8x10 'S. 


670.86341454 

l.OxlO' 5 ?. 


670.86345124 

9.1x10' 5 % 


I 


5 


1002.15520225 


1002.19978814 


1002.18735784 


1002.19300449 


\ 

l 

1 




j 


4'. 7x10 -3 % 


3.2xl0'\ 


3.8x10" 3 % 


ll 

r 

< 

\ 

i 

i 


6 


1399.70436269 


i 


1399.59981662 
7 . 5x10 _3 % 


1399.64448571 

4.3x10 -3 % 


) 

i 


7 

| 


1863.51172599 


T " ' r ■ ■ 1 

i 




1959.59556134 
5.1 % 



Mode 

No. 


Freq. (rad. sec) and difference in % 


Delta m. (S) 


P-M m.-A(S) 


M-D m.(D)*** 


M-D m.(S)*** 


1 . 


75.09994535 
0 % 


75.09994865 

£ 

4.0x10 % 


75.09994539 
0 7. 


75.09994534 
0 7. 


2 


207.01589228 207.01589071 

5.0x10' 8 % 0 7 , 


207.01588836 
1.5xl0' 6 7. 


207.01588891 
0 7. 


3 


405.83398145 405.83403455 

1.5x10' 6 % 2. lxlO -5 % 


405.83395151 
9.5xl0' 6 7. 


405.83401120 
2. 0x10' 5 7. 


4 


670.86345124 

7.5x10 _5 % 


670.86256288 
2.2xl0 -4 % 


670.86374475 
5. lxlO -5 7. 


670.86118054 
4.3xl0 -4 7. 


5 




1002.23499946 
8.0x10" 3 7 , 


1002.19909306 
4. 4x10 " 4 7. 


1002.21779856 
6. 3xl0 3 % 


6 




1399.17084415 
4. 5x10" 2 7. 


1399.94044241 
3. 4x10' 2 7. 





★Sources of comparison frequencies for system 1 through 3 are given in E-2-2. 
(Additional footnotes on next page). 
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**In all the tables which follow, simple and evident abbreviations are used. 
For example, P-M m.-A(D) means the P-M method, sub-procedure A in double 
precision arithmetic. 

***Modified System 1. 
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b. System 2. (in-plane) 



Mode . 
NO 


Freq. ( rad/sec) and difference in % 


- 


Comparison values 


Delta m. (d) 


P-M m.-A (D) 


P-M m.-P (D) 

11.80213746 
1.4xl0~ 5 % 


1 


11.80213386 


11.80213389 

5. 5x10" 5 % 


11.80214614 
9.5xlO" 5 $> 


2 


73.96272297 


73.96 272167 
1.7xl0 -6 % 


73.96272478 
1.7x10-6 £ 


73.96272367 
0 % 


3 


2 07.09776595 


207.09776685 

4.1xl0“ 8 $ 


207.09776388 
l.lxlO" 6 # 


207.09776562 
0 # 




4 


405.82896584 


405.82894132 
6. 1x10 " 6 # 


405.82893748 
7.0xlO~ 6 % 


405.82900191 
9.1xl0* 6 % 


5 • • 


- 670.86435913 


670.86200817 
2.2x10-5 % 


670.86326697 
1.7xlO~ 5 % 


670.86325735 
1.7x10-5 % 




6 


1002.1551877 


1002.26833752 

1.1x10-2 % 


1002.18809484 
4.4xl0“3 $ 


1002.19715557 
4.2xl0" 3 # . 


7 


1399.70436329 


1399.08849534 
•2.0xl0 -2 # 


1399.61231065 

6. 6x10“ 3 % 


1399.14599562 
4.2x10-3 % 


8 


1863.51172626 




1862.55983728 
5.1xlO" 2 % 


i 

1 

1 



c. System 3. (in-plane) 



Mode 

NO 


Freq. ( rad/sec) and difference in % 


Comparison values 


Delta m. (d) 


P-M m.-A (D) 


P-M m.-P (D) 


1 


51.75397285 


51.75397263 
0 % 


51.75397569 
0 $> 


51.75397564 
0 % 


2 


167.71502090 


167.71602009 
0 % 


167.71602200 

0 


167.71602189 
0 % 


3 


349.92609071 


349.92608231 
l.lxlO" 7 % 


349.92610259 
3.4xl0“ 6 % 


349.92610209 

3.4x10-6 % 


4 


598.39432089 


598.39381063 
8.5x10-5 % 


598.39446341 
2.5xlO- 5 % 


598.39407344 
4. 2x10" 5 # 


5 


913.12074582 


913.13150640 
1. 2xlO~ 3 % 


913.13180993 
1.2x10" 3 % 


913.13109335 
1.2xlO- 3 % 


6 


1294.10536551 


1293.77838653 
2.7xlO” 2 % 


1294.11223844 
5.8xl0" 4 % 


1294.11761177 
1.9xlO" 4 # 


7 


1741^34817946 




1748.59011823 
4. 2xl0 -1 % 
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E-2-2 Sources of comparison frequencies for system 1 through 3. 

The governing fourth order equation of a homogeneous single component 
straight system 

y v _ ajL-.'i 

h EX J' 

was solved to give 

y ■ a cos d x + b cosheX x + c sin ^ sintvXx 

Upon applying boundary conditions for System 1 (y ■ y' =0 at x * 0, 
x « L) and simplifying, it was found that the eigenvalues of the system 
must satisfy 

cos o/- cosh(X “ 1 (E-2-2-1 ) 

where _ Vx tQ 4 ^ 

^ J EX (E-2-2-2) 

Upon applying the boundary conditions for System 2(y = y' * 0 at x * 0, 
y" = y'" = o at x ■ L), it was found that the eigenvalues of the system 
must satisfy 

cos U * cosho(" -1 (E - 2- 2- 3) 

Upon applying the boundary conditions of System 3(y • y' ■ 0 at x » 0, 
y a y" = 0 at x = L)> it was found that the eigenvalues must satisfy 
sino( • coshO(- cos o( • cosh oi * 0 (E-2-2-4) 

The roots of equation E-2-2-1, E-2-2-3 and E-2-2-4, from which the natural 
frequencies of the corresponding systems may easily be found, were ob- 
tained through iteration using a digital computer. 
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E-2-3 Curved section (System 4) 

System parameters: 




End conditions: 



radius of curvature 
included angle of arc 
diameter 
wall thickness 



50 inches 
270 degrees 
2.0 inches 
.125 inches 



Intensive properties: 



density 
shear modulus 
elastic modulus 



556. lbs/cu-ft 
12 x 106 psi 
30xl0 6 psi 



fixed-fixed. 



Mathematical model: 



lumped mass with shear deflection and rota- 
tional inertia neglected. 



Results. 



a. In-plane vibration* 



Mode 


Frequency (rad/sec) 


! No. 


Delta m. (D) 


M-D m. (D) 


P-M m.-A (D) 


P-M m.-C (D) 


1 


66.69934902 


66.69934902 


66.69934980 


66.69934948 


2 


171.22352538 


171.22352538 


171.22352610 


171.22352592 


i 3 

» 


335.35007083 


335.35007083 




335.35007083 


335.35007174 


] 4 




535.64529414 


535.6429665 


535.64529637 


5 




773.44358487 


773.44358842 


773.44358812 



b. Out-of-plane vibration* 



Mode 

No. 


Frequency (rad/sec) 


Delta m. (D) 


M-D m. (D) 


P-M m.-A (D) 


1 


35.36309735 


35.36309735 


35.36309778 


2 


94.42588921 


94.42588921 


94.42588653 


3 


207.89287257 


207.89278647 


207.89287345 I 

. ...... j 


4 


369.11199342 


369.11278647 


369.11278689 j 

j 


5 




571.22235666 


571.22235577 



*See footnote next page 
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*The fundamental frequency value of in-plane and out-of-plane case are ^ 4 . 
close to the values obtained by employing equation O=‘p(o 0 *( E^/CuR ) 
after obtaining F( o< ) from the curves in reference [5]. Natural 
frequency values from the three methods are the same up to six digits. 

One can get more higher mode frequencies using the M-D and P-M method. 
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E-3 Other Sample Problems 
A, System 5 



System parameters 



Parameter 


component 


i 


2 


3 


4 


5 


6 


7 


diameter (in) 


2.0 


2.0 


2.0 


H 


2.0 


2.0 


2.0 










A 








wall thickness 


). 75 


0.75 


0.75 


N 


0.75 


0.75 


0.75 


(in) 








G 








length (in) 


130. 


1.40 


6.0 


E 


60.21 


1.50 


6.28 










R 








radius of curve 




4.0 






30. 




4.0 


(in) 
















incl. angle of arc 




20. 






115. 




90. 



Intensive properties of all components except hangers: 

density: 556 lbs/cu-ft.* 

shear modulus: 12x10^ psi 
elastic modulus: 30x10^ psi 



Hanger spring rates: 

CLX , CLY , CLZ 1000 lbs/ in. 

CTX , CTY , CTZ 1000 in- lb/rad. 

End conditions: 



both ends of system fixed. 

Mathematical model: Distributed mass with the effects of shear deflec- 

tion and rotational inertia considered. 



t 



1 








*This fictitious value was chosen to permit comparing with results of 
Fink [4] who used the same value. 
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Results, (in-plane vibration) 



Mode 

NO 




Frequency ( rad/ sec ) 




Delta m. (d) 


M-D o. (D) 


P-M m.-A (D) 


P-M m.-C (D) 


1 


116.89479180 








2 


304.88354398 








3 

i 


566.43670976 








4 




826.8828 2292 


826.88281950 


826.88281929 


5 

L 




1126.83773518 


1126.83799669 


1126.83797508 


6 




1666.58935151 


1666.58812881 


1666.58878192 


7 




2099.54173011 


2099.5352509 


2099.52844030 



Results in single precision arithmetic and comparison with double precision results. 
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B. System 6 (single branch system) 



System parameters: 





components 


parameter 


i 


2 


3 


4 


dia. (inch) 


2.0 


2.0 


2.0 


2.0 


wall thickness 


.125 


.125 


.125 


.125 


(inch) 

length (inch) 


100. 


100. 


.01 


100. 


incl. angle of 
arc. (degree) 


0 


0 


45 


0 




Intensive properties: 

density 556 lbs/cu-ft 

shear modulus 12x10** psi 

elastic 30x10^ psi 

modulus 

End conditions: All ends fixed. 

Mathematical model: Distributed mass with the effects of shear deflec- 

tion and rotational inertia neglected. 

Results 



Mode 

No. 


Frequency (rad/ sec) 




Delta m. (D) M-D m. (D) 


P-M m.-A (D) 


P-M m.-C (D) 


1 


194.25470558 




t 

i 

i 


2 


278.61078241 






3 


628.59556355 






4 


747.68138190 747.68087208 


747.68137310 


747.68137759 


5 


1300.29478437 


1300.31646085 


1300.31641281 


6 


1381.79339427 


1381.86234251 


1381.86230734 


7 


1991.80356348 


1991.38946110 


1991.38946101 


8 


2257.92926800 

1 — J 


2257.58161443 

— i 


2257.58016324 
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